Coverage for matrix.py : 63%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
#! /usr/bin/env python # -*- coding: utf-8 -*-
This file defines a class :class:`Matrix`, designed to be as complete as possible. *Do not worry, I was not asking you to do as much.*
Examples --------
Importing the module:
>>> from matrix import * >>> from matrix import Matrix as M # shortcut
Defining a matrix by giving its list of rows:
>>> A = M([[1, 0], [0, 1]]) >>> A == eye(A.n) True >>> B = 2*(A**2) + 4*A + eye(A.n) >>> B [[7, 0], [0, 7]] >>> B == 7 * eye(A.n) True
Indexing and slicing:
>>> A[1,:] = 2; A [[1, 0], [2, 2]] >>> A[0, 0] = -5; A [[-5, 0], [2, 2]]
Addition, multiplication, power etc:
>>> C = eye(2); C [[1, 0], [0, 1]] >>> C + (3 * C) - C [[3, 0], [0, 3]] >>> (4 * C) ** 2 [[16, 0], [0, 16]]
Many more examples are given below:
-----------------------------------------------------------------------------
Things that could still be worked on for this solution ------------------------------------------------------ .. todo:: Implement the **QR**, **SVD** and other **matrix decompositions**. .. todo:: Try to add a randomized matrix decomposition (or any *less-original* matrix decomposition method)? Note: I worked on this aspect, for a project in January 2016 for my M.Sc. : `<https://bitbucket.org/lbesson/mva15-project-parcimonie-compressed-sensing/>`_. .. todo:: Implement a nice wrapper for a linear equations solver (with LU). .. todo:: More doctests for :py:func:`PLUdecomposition`, and implement the non-permuted LU decomposition? .. todo:: Add more doctests and examples for Gauss, Gauss-Jordan, Gram-Schmidt (:py:func:`gauss`, :py:func:`gauss_jordan`, :py:func:`gram_schmidt`)?
.. note:: Interactive examples?
See the other file `tests.py <tests.html>`_ for *many* examples.
- *Date:* Saturday 18 juin 2016, 10:31:25. - *Author:* `Lilian Besson <https://bitbucket.org/lbesson/>`_ for the `CS101 course <http://perso.crans.org/besson/cs101/>`_ at `Mahindra Ecole Centrale <http://www.mahindraecolecentrale.edu.in/>`_, 2015, - *Licence:* `MIT Licence <http://lbesson.mit-license.org>`_.
.. seealso::
I also wrote a complete solution for the other project I was in charge of, `about numerical algorithms to compute integrals <http://mec-cs101-integrals.readthedocs.io/en/latest/integrals.html>`_. """
# Experimental: computation with decimal numbers to improve decimal precision
""" Extended :class:`decimal.Decimal` class to improve the ``str`` and ``repr`` methods.
If there is not digit after the comma, print it as an integer. """ if int(self) == self: return "{}".format(int(self)) else: return _Decimal.__str__(self)
# Experimental: computation with fraction to be exact and not numerically approximative !
""" Extended :class:`fractions.Fraction` class to improve the ``str`` and ``repr`` methods.
If the denominator is 1, print it as an integer. """ else:
# ======================================================================== """ A class to represent matrices of size ``(n, m)``.
``M = Matrix(listrows)`` will have three attributes:
- :py:data:`M.listrows` list of rows vectors (as list), - :py:data:`M.n` or :py:data:`M.rows` number of rows, - :py:data:`M.` or :py:data:`M.cols` number of columns (ie. length of the rows).
All the required special methods are implemented, so :class:`Matrix` objects can be used as numbers, with a very natural syntax.
.. warning:: All the rows should have the same size. """
""" Create a :class:`Matrix` object from the list of row vectors ``M``.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> A.listrows [[1, 2, 3], [4, 5, 6]] """ self._m = 0 else: # Now we do get a fresh copy of that list. #: self.listrows is the list of rows for self # FIXME We should forbid modifying these attributes from outside the class except: raise ValueError("Matrix() accepts only a list of rows vectors (ie. list of lists) as its argument.")
# This decorator @property makes this method an attributes # cf. https://docs.python.org/2/library/functions.html#property def n(self): """ Getter for the read-only attribute ``A.n`` (number of rows).
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> A.n 2 >>> A.rows == A.n True """
# This decorator @property makes this method an attributes # cf. https://docs.python.org/2/library/functions.html#property def m(self): """ Getter for the read-only attribute ``A.m`` (size of the rows, ie. number of columns).
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> A.m 3 >>> A.cols == A.m True """
# ==================================================================== # Methods for reading and accessing
""" ``A[i, j]`` <-> ``A.listrows[i][j]`` reads the (``i, j``) element of the matrix ``A``.
- *Experimental* support of slices: ``A[a:b:k, j]``, or ``A[i, c:d:l]`` or ``A[a:b:k, c:d:l]``. - Default values for ``a`` and ``c`` is a **start point** of ``0``, ``b`` and ``d`` is a **end point** of maximum size, and ``k`` and ``l`` is a **step** of ``1``.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> A[0, 0] 1 >>> A[0, :] [[1, 2, 3]] >>> A[-1, :] [[4, 5, 6]] >>> A[:, 0] [[1], [4]] >>> A[1:, 1:] [[5, 6]] >>> A[:, ::2] [[1, 3], [4, 6]] """ # i is an integer, j is a slice object # i is a slice object, j is an integer # i and j are a slice objects # In case i and j are neither integers nor slice objects raise ValueError("Matrix.__getitem__ invalid argument. A[i, j] with i = {} (type(i) is {}) and j = {} (type(i) is {}).".format(i, type(i), j, type(j)))
""" ``A[i, j] = value``: will update the ``(i, j)`` element of the matrix ``A``.
- Support for slice arguments: ``A[a:b:k, j] = sub_row``, or ``A[i, c:d:l] = sub_column`` or ``A[a:b:k, c:d:l] = submatrix``. - Default values for ``a`` and ``c`` is a **start point** of ``0``, ``b`` and ``d`` is a **end point** of maximum size, and ``k`` and ``l`` is a **step** of ``1``.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> A[0, 0] = 4; A [[4, 2, 3], [4, 5, 6]] >>> A[:, 0] [[4], [4]] >>> A[-1, :] = 9; A [[4, 2, 3], [9, 9, 9]] >>> A[1, 1] = 3; A [[4, 2, 3], [9, 3, 9]] >>> A[0, :] = [3, 2, 1]; A [[3, 2, 1], [9, 3, 9]] >>> A[1:, 1:] = -1; A [[3, 2, 1], [9, -1, -1]] >>> A[1:, 1:] *= -8; A [[3, 2, 1], [9, 8, 8]] """ # This is the simple case, the one we use the most # i is an integer, j is a slice object # for j0 in _slice_to_range(j): except Exception: try: self.listrows[i][j0] = value[j_value] # list except Exception: self.listrows[i][j0] = value # End for loop j0 else: fail = True # i is a slice object, j is an integer # for i0 in _slice_to_range(i): # End for loop i0 elif isinstance(j, slice): # i and j are a slice objects i_value = 0 j_value = 0 # for i0 in _slice_to_range(i): for i0 in range(*i.indices(self.rows)): # for j0 in _slice_to_range(j): for j0 in range(*j.indices(self.cols)): try: self.listrows[i0][j0] = value[i_value][j_value] # sub-matrix except Exception: try: self.listrows[i0][j0] = value[i_value] # list except Exception: self.listrows[i0][j0] = value j_value += 1 # End for loop i0 i_value += 1 # End for loop i0 else: fail = True # In case i and j are neither integers nor slice objects raise ValueError("Matrix.__setitem__ invalid argument. A[i, j] with i = {} (type(i) is {}) and j = {} (type(i) is {}).".format(i, type(i), j, type(j)))
# row and col to access a row or a column
""" ``A.row(i)`` <-> *extracts* the ``i``-th row of ``A``, as a *new* matrix.
.. warning:: Modifying ``A.row(i)`` does NOT modify the matrix ``A``.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> A.row(0) [[1, 2, 3]] >>> A.row(1) [[4, 5, 6]] >>> r = A.row(0); r *= 3 >>> A # it has not been modified! [[1, 2, 3], [4, 5, 6]] """
""" ``A.col(j)`` <-> *extracts* the ``j``-th column of ``A``, as a new matrix.
.. warning:: Modifying ``A.col(j)`` does NOT modify the matrix A.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> A.col(0) [[1], [4]] >>> A.col(2) [[3], [6]] >>> c = A.col(1); c *= 6 >>> A # it has not been modified! [[1, 2, 3], [4, 5, 6]] """
# ==================================================================== # Special method for copying (not required in the project)
#def __hash__(self): # """ hash(A) <-> A.__hash__() computes the hash of the matrix (just depends on A.listrows). # # - Is required if we want to be able to insert a matrix in a set or dictionary. # - FIXME: unhashable type: 'list' # """ # return hash(self.listrows)
""" ``A.copy()`` <-> a shallow copy of the matrix ``A`` (ie. a new and fresh matrix with same values).
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> B = A.copy() >>> A[0, 0] = -10; A [[-10, 2, 3], [4, 5, 6]] >>> B # It has not been modified! [[1, 2, 3], [4, 5, 6]] """
# ==================================================================== # Length and shape
""" ``len(A)`` returns ``A.n * A.m``, the number of values in the matrix.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> len(A) 6 >>> len(A) == A.n * A.m True """
# This decorator @property makes this method an attributes # cf. https://docs.python.org/2/library/functions.html#property def shape(self): """ ``A.shape`` is ``(A.n, A.m)`` (similar to the shape attribute of NumPy arrays).
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> A.shape (2, 3) """
# ==================================================================== # Transposition
""" ``A.transpose()`` is the transposition of the matrix ``A``.
- Returns a new matrix! - Definition: if ``B = A.transpose()``, then ``B[i, j] is A[j, i]``.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> A.transpose() [[1, 4], [2, 5], [3, 6]] >>> A.transpose().transpose() == A True """
def T(self): """ ``A.T`` <-> ``A.transpose()`` is the transposition of the matrix ``A``, useful shortcut as in NumPy.
>>> B = Matrix([[1, 4], [2, 5], [3, 6]]) >>> B.T [[1, 2, 3], [4, 5, 6]] >>> B == B.T.T True """
# ==================================================================== # Methods for pretty-printing
""" ``str(A)`` <-> ``A.__str__()`` converts the matrix ``A`` to a string (showing the list of rows vectors).
>>> B = Matrix([[1, 4], [2, 5], [3, 6]]) >>> str(B) '[[1, 4], [2, 5], [3, 6]]' """ except Exception: str(self.listrows)
""" ``repr(A)`` <-> ``A.__repr__()`` converts the matrix A to a string (showing the list of rows vectors).
>>> B = Matrix([[1, 4], [2, 5], [3, 6]]) >>> repr(B) '[[1, 4], [2, 5], [3, 6]]' """ return str(self)
# Comparing ==
r""" ``A == B`` <-> ``A.__eq__(B)`` compares the matrix ``A`` with ``B``.
- Time complexity is :math:`\mathcal{O}(n m)` for matrices of size ``(n, m)``.
>>> B = Matrix([[1, 4], [2, 5], [3, 6]]) >>> B == B True >>> B + B + B == 3*B == B + 2*B == 2*B + B True >>> B - B + B == 1*B == -B + 2*B == 2*B - B == 2*B + (-B) True >>> B != B False """ # return all(a == b for a, b in zip(self, B)) else: return False except Exception: return False
r""" ``A.almosteq(B)`` compares the matrix ``A`` with ``B``, numerically with an error threshold of ``epsilon``.
- Default epsilon is :math:`10^{-10}`. - Time complexity is :math:`\mathcal{O}(n m)` for matrices of size ``(n, m)``.
>>> B = Matrix([[1, 4], [2, 5], [3, 6]]) >>> C = B.copy(); C[0,0] += 4*1e-6 >>> B == C False >>> B.almosteq(C) False >>> B.almosteq(C, epsilon=1e-4) True >>> B.almosteq(C, epsilon=1e-5) True >>> B.almosteq(C, epsilon=1e-6) False """ else: return False except Exception: return False
# Comparing <
r""" ``A < B`` <-> :math:`A_{i,j} < B_{i,j} \forall i,j` compares the matrix ``A`` with ``B``.
- Time complexity is :math:`\mathcal{O}(n m)` for matrices of size ``(n, m)``. - Time complexity is :math:`\mathcal{O}(n m)` for matrices of size ``(n, m)``. - ``A > B``, ``A <= B``, ``A >= B`` are all computed automatically with :py:meth:`__eq__` and :py:meth:`__lt__`.
>>> B = Matrix([[1, 4], [2, 5], [3, 6]]) >>> B < B False >>> B < B + 4 True >>> B > B False >>> B > B - 12 True """ try: if self.n == B.n and self.m == B.m: return all(self[i, j] < B[i, j] for j in range(self.m) for i in range(self.n)) # return all(a < b for a, b in zip(self, B)) else: return False except Exception: return False
# ==================================================================== # Methods for computing
# Sum (left and right)
r""" ``A + B`` <-> ``A.__add__(B)`` computes the sum of the matrix ``A`` and ``B``.
- Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m)` for matrices of size ``(n, m)``. - If ``B`` is a number, the sum is done coefficient wise.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> A + A [[2, 4, 6], [8, 10, 12]] >>> B = ones(A.n, A.m); B [[1, 1, 1], [1, 1, 1]] >>> A + B [[2, 3, 4], [5, 6, 7]] >>> B + A [[2, 3, 4], [5, 6, 7]] >>> B + B + B + B + B + B + B [[7, 7, 7], [7, 7, 7]] >>> B + 4 # Coefficient wise! [[5, 5, 5], [5, 5, 5]] >>> B + (-2) # Coefficient wise! [[-1, -1, -1], [-1, -1, -1]] >>> B + (-1.0) # Coefficient wise! [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] """ # Sum of two matrices else: # Sum of matrix A and a number B
r""" ``B + A`` <-> ``A.__radd__(B)`` computes the sum of ``B`` and the matrix ``A``.
- Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m)` for matrices of size ``(n, m)``. - If ``B`` is a number, the sum is done coefficient wise.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> 1 + A [[2, 3, 4], [5, 6, 7]] >>> B = ones(A.n, A.m) >>> 4 + B # Coefficient wise! [[5, 5, 5], [5, 5, 5]] >>> (-2) + B # Coefficient wise! [[-1, -1, -1], [-1, -1, -1]] >>> (-1.0) + B # Coefficient wise! [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] """ # Sum of two matrices # (never used here : B + A <-> B.__add__(A)) assert self.n == B.n and self.m == B.m return Matrix([[B[i, j] + self[i, j] for j in range(self.m)] for i in range(self.n)]) else: # Sum of matrix A and a number B (coefficients wise)
# ==================================================================== # Substraction (left and right)
r""" ``A - B`` <-> ``A.__sub__(B)`` computes the difference of the matrix ``A`` and ``B``.
- Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m)` for matrices of size ``(n, m)``. - If ``B`` is a number, the sum is done coefficient wise.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> B = ones(A.n, A.m) >>> A - B [[0, 1, 2], [3, 4, 5]] >>> B - A [[0, -1, -2], [-3, -4, -5]] >>> A - 1 # Coefficient wise! [[0, 1, 2], [3, 4, 5]] >>> B - 2 # Coefficient wise! [[-1, -1, -1], [-1, -1, -1]] >>> (A - 3.14).round() # Coefficient wise! [[-2.14, -1.14, -0.14], [0.86, 1.86, 2.86]] """ # Sum of two matrices else: # Sum of matrix A and a number B
r""" ``-A`` <-> ``A.__neg__()`` computes the opposite of the matrix ``A``.
- Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m)` for a matrix of size ``(n, m)``.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> -A [[-1, -2, -3], [-4, -5, -6]] >>> A - A == A + (-A) True >>> -(-A) == A True >>> -------A == -A # Crazy syntax! True >>> s = '-------' >>> len(s) % 2 == 1 # We check that we had an od number of minus symbol True """
r""" ``+`` <-> ``A.__pos__()`` computes the positive of the matrix A.
- Returns a new matrix! - Useless? - Time and memory complexity is :math:`\mathcal{O}(n m)` for a matrix of size ``(n, m)``.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> +A == A True >>> +-+-+-+-+++----+-+-+----++++A == A # Crazy syntax, again! True >>> s = '+-+-+-+-+++----+-+-+----++++' >>> s.count('-') % 2 == 0 # We check that we had an even number of minus symbol True """
r""" ``B - A`` <-> ``A.__rsub__(B)`` computes the difference of ``B`` and the matrix ``A``.
- Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m)` for matrices of size ``(n, m)``. - If ``B`` is a number, the sum is done coefficient wise. - If ``B`` is a :class:`Matrix` object, ``B - A`` will in fact be ``B.__sub__(A)`` and not ``A.__rsub__(B)``.
>>> A = Matrix([[1, 2, 3], [4, 5, 6]]) >>> 1 - A # Coefficient wise! [[0, -1, -2], [-3, -4, -5]] >>> B = ones(A.n, A.m) >>> (-1) - B # Coefficient wise! [[-2, -2, -2], [-2, -2, -2]] >>> ((-1) - B) == -(1 + B) == -(B + B) True """ # Sum of two matrices # (never used here : B - A <-> B.__sub__(A)) assert self.n == B.n and self.m == B.m return Matrix([[B[i, j] - self[i, j] for j in range(self.m)] for i in range(self.n)]) else: # Sum of matrix A and a number B (coefficients wise)
# ==================================================================== # Product (left and right)
r""" ``A * B`` <-> ``A.__mul__(B)`` computes the product of the matrix ``A`` and ``B``.
- Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m p)` for a matrix ``A`` of size ``(n, m)`` and ``B`` of size ``(m, p)``. - If ``B`` is a number, the product is done coefficient wise.
.. warning:: Matrix product is not commutative!
>>> A = Matrix([[1, 2], [3, 4]]) >>> B = eye(A.n); B [[1, 0], [0, 1]] >>> A * B == B * A == A True >>> A * A [[7, 10], [15, 22]] >>> A * (A * A) == (A * A) * A True >>> A * 1 == A # Coefficient wise! True >>> A * 12.011993 # Coefficient wise! [[12.011993, 24.023986], [36.035979, 48.047972]] """ # Product of two matrices else: # Product of matrix A and a number B (coefficients wise)
r""" ``B * A`` <-> ``A.__rmul__(B)`` computes the product of ``B`` and the matrix ``A``.
- Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m p)` for a matrix ``A`` of size ``(n, m)`` and ``B`` of size ``(m, p)``. - If B is a number, the product is done coefficient wise. - If ``B`` is a :class:`Matrix` object, ``B * A`` will in fact be ``B.__mul__(A)`` and not ``A.__rmul__(B)``.
.. warning:: Matrix product is not commutative!
>>> A = Matrix([[1, 2], [3, 4]]) >>> 1 * A == A # Coefficient wise! True >>> 12.011993 * A # Coefficient wise! [[12.011993, 24.023986], [36.035979, 48.047972]] """ # Product of two matrices # (never used here : B * A <-> B.__mul__(A)) assert self.n == B.m return Matrix([[sum(B[i, k] * self[k, j] for k in range(B.m)) for j in range(self.m)] for i in range(B.n)]) else: # Product of matrix A and a number B (coefficients wise)
r""" ``A.multiply_elementwise(B)`` computes the product of the matrix ``A`` and ``B``, element-wise (it is called a **Hadamard product**).
- Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m p)` for a matrix ``A`` of size ``(n, m)`` and ``B`` of size ``(m, p)``.
>>> A = Matrix([[1, 2], [3, 4]]) >>> B = eye(A.n) >>> A.multiply_elementwise(B) [[1, 0], [0, 4]] >>> A.multiply_elementwise(A) # A .^ 2 in Matlab? [[1, 4], [9, 16]] """ raise ValueError("A.multiply_elementwise(B): B has to be a Matrix object.") else:
# ==================================================================== # Division (left and right)
r""" ``A / B`` <-> ``A * (B ** (-1))`` computes the division of the matrix ``A`` by ``B``.
- Returns a new matrix! - Performs **true division**! - Time and memory complexity is :math:`\mathcal{O}(n m p \max(m, p)^2)` for a matrix ``A`` of size ``(n, m)`` and ``B`` of size ``(m, p)``. - If ``B`` is a number, the division is done coefficient wise.
>>> A = Matrix([[1, 2], [3, 4]]) >>> B = eye(A.n) >>> B.almosteq(A / A) True >>> C = B.map(float) >>> A / C == A * C == A True >>> A / B == A * B == A True >>> A / 2 # Coefficient wise! [[0.5, 1.0], [1.5, 2.0]] >>> A / 2.0 # Coefficient wise! [[0.5, 1.0], [1.5, 2.0]] """ # print("self.__div__:", B, type(B)) # DEBUG. # print("self.__div__:", B, type(B)) # DEBUG. # Product of two matrices return self * (B.inv()) else: # Division of matrix A and a number B (coefficients wise) # return Matrix([[self[i, j] / float(B) for j in range(self.m)] for i in range(self.n)])
r""" ``A // B`` <-> ``A * (B ** (-1))`` computes the division of the matrix ``A`` by ``B``.
- Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m p)` for a matrix ``A`` of size ``(n, m)`` and ``B`` of size ``(m, p)``. - If ``B`` is a number, the division is done coefficient wise with an **integer division** ``//``.
>>> A = Matrix([[1, 2], [3, 4]]) >>> B = eye(A.n); C = B.map(float) >>> A // C == A * C == A True >>> A // B == A * B == A True >>> A // 2 # Coefficient wise! [[0, 1], [1, 2]] >>> A // 2.0 # Coefficient wise! [[0.0, 1.0], [1.0, 2.0]] """ if isinstance(B, Matrix): # Product of two matrices return self * (B.inv()) else: # Division of matrix A and a number B (coefficients wise) return Matrix([[self[i, j] // B for j in range(self.m)] for i in range(self.n)])
r""" ``A % b`` <-> ``A.__mod__(b)`` computes the modulus coefficient-wise of the matrix ``A`` by ``b``.
- Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m)` for a matrix ``A`` of size ``(n, m)``.
>>> A = Matrix([[1, 2], [3, 4]]) >>> A % 2 [[1, 0], [1, 0]] >>> (A*100) % 31 [[7, 14], [21, 28]] >>> (A*100) % 33 == A # Curious property True >>> (A*100) % 35 [[30, 25], [20, 15]]
.. warning:: ``A % B`` for two matrices means the coefficient-wise modulus.
>>> A = Matrix([[1, 2], [3, 4]]) >>> B = Matrix([[2, 3], [2, 2]]) >>> A % B [[1, 2], [1, 0]] """ # Product of two matrices return Matrix([[self[i, j] % b[i, j] for j in range(self.m)] for i in range(self.n)]) else:
r""" ``B / A`` <-> ``A.__rdiv__(B)`` computes the division of ``B`` by ``A``.
.. warning:: If ``B`` is ``1`` (``B == 1``), ``1 / A`` is ``A.inv()`` (special case!)
- If ``B`` is a number, the division is done coefficient wise. - Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m p)` for a matrix ``A`` of size ``(n, m)`` and ``B`` of size ``(m, p)``.
>>> A = Matrix([[1, 2], [3, 4]]) >>> Ainv = Matrix([[-2.0, 1.0], [1.5, -0.5]]) >>> B = eye(A.n) >>> B == A * Ainv == Ainv * A True >>> 1 / B == B == B / 1 True >>> C = B.map(float) >>> 1 / B == B == B / 1 True >>> A.inv() == 1 / A # special case! True >>> 1 / A # This is like 1 / A [[-2.0, 1.0], [1.5, -0.5]] >>> 2 / (2*A) # Warning This is coefficient wise ! # doctest: +ELLIPSIS [[1.0, 0.5], [0.333333..., 0.25]] """ # print("self.__rdiv__:", B, type(B)) # DEBUG. if B == 1: return self.inv() elif isinstance(B, Matrix): return B * (self.inv()) else: # Division of a number B and matrix A (coefficients wise) return Matrix([[B / self[i, j] for j in range(self.m)] for i in range(self.n)])
r""" ``B // A`` <-> ``A.__rdiv__(B)`` computes the division of ``B`` by ``A``.
.. warning:: If ``B`` is ``1`` (``B == 1``), ``1 / A`` is ``A.inv()`` (special case!)
- If ``B`` is a number, the division is done coefficient wise. - Returns a new matrix! - Time and memory complexity is :math:`\mathcal{O}(n m p)` for a matrix ``A`` of size ``(n, m)`` and ``B`` of size ``(m, p)``.
>>> A = Matrix([[1, 2], [3, 4]]) >>> B = eye(A.n) >>> 1 // B == B == B // 1 True >>> C = B.map(float) >>> 1 // B == B == B // 1 True >>> A.inv() == 1 // A # special case! True >>> 2 // (2*A) # XXX This is coefficient wise ! [[1, 0], [0, 0]] """ # print("self.__rdiv__:", B, type(B)) # DEBUG. if B == 1: return self.inv() elif isinstance(B, Matrix): return B * (self.inv()) else: # Division of a number B and matrix A (coefficients wise) return Matrix([[B // self[i, j] for j in range(self.m)] for i in range(self.n)])
# ==================================================================== # Power, exponential and inverse
r""" ``A ** k`` <-> ``A.__pow__(k)`` to compute the product of the square matrix ``A`` (with the quick exponentation trick).
- Returns a new matrix! - ``k`` has to be an integer (``ValueError`` will be returned otherwise). - Time complexity is :math:`\mathcal{O}(n^3 \log(k))` for a matrix ``A`` of size (n, n). - Memory complexity is :math:`\mathcal{O}(n^2)`. - It uses ``A.inv()`` (:py:meth:`inv`) to (try to) compute the inverse if ``k < 0``. - More details are in `the solution for the Problem II of the 2nd Mid-Term Exam for CS101 <http://perso.crans.org/besson/cs101/Exams/Second_MidTerm_Exam/>`_.
>>> A = Matrix([[1, 2], [3, 4]]) >>> A ** 1 == A True >>> A ** 2 [[7, 10], [15, 22]] >>> A * A == A ** 2 True >>> B = eye(A.n) >>> B == B ** 1 == A ** 0 == B ** 0 True >>> divmod(2015, 2) (1007, 1) >>> 2015 == 1007*2 + 1 True >>> A ** 2015 == ((A ** 1007) ** 2 ) * A True >>> C = diag([1, 4]) >>> C ** 100 [[1, 0], [0, 1606938044258990275541962092341162602522202993782792835301376]] >>> C ** 100 == diag([1**100, 4**100]) True
It also accept negative integers:
>>> A ** (-1) == A.inv() True >>> C = (A ** (-1)); C [[-2.0, 1.0], [1.5, -0.5]] >>> C * A == eye(A.n) == A * C True >>> C.listrows # Rounding mistakes can happen (but not here) [[-2.0, 1.0], [1.5, -0.5]] >>> D = C.round(); D.listrows [[-2.0, 1.0], [1.5, -0.5]] >>> D * A == eye(A.n) == A * D # No rounding mistake! True >>> (C * A).almosteq(eye(A.n)) True >>> (A ** (-5)) == (A ** 5).inv() == (A.inv()) ** 5 False >>> (A ** (-5)).round() == ((A ** 5).inv()).round() == ((A.inv()) ** 5).round() # No rounding mistake! True """ raise ValueError("A ** k: k should be an integer (here k = {}).".format(k)) # A ^ k = (A ^ (-1)) ^ (-k) return (self.inv()) ** (abs(k)) return eye(self.n) # This is a convention : a ** 0 = 1_n,n # Remark: this case is not tail recursive (we could have used an accumulator) else: # XXX should never happen! raise ValueError("A ** k: invalid value for the power k = {k}.".format(k=k))
r""" ``A.exp()`` computes a numerical approximation of the exponential of the square matrix ``A``.
- Raise a ValueError exception if ``A`` is not square. - Note: :math:`\exp(A) = \mathrm{e}^A` is defined as the series :math:`\sum\limits_{k=0}^{+\infty} \frac{A^k}{k!}`. - We only compute the first ``limit`` terms of this series, hopping that the partial sum will be close to the entire series. - Default value for ``limit`` is 30 (it should be enough for any matrix).
>>> import math >>> e = math.e >>> I = eye(10); I[0, :] [[1, 0, 0, 0, 0, 0, 0, 0, 0, 0]] >>> I * e == I.exp() == diag([e] * I.n) # Rounding mistakes! False >>> (I * e).round() == I.exp().round() == diag([e] * I.n).round() # No more rounding mistakes! True >>> C = diag([1, 4]) >>> C.exp() == diag([e ** 1, e ** 4]) == diag([math.exp(1), math.exp(4)]) # Rounding mistakes! False >>> C.exp().almosteq(diag([e ** 1, e ** 4])) # No more rounding mistakes! True >>> diag([e ** 1, e ** 4]).almosteq(diag([math.exp(1), math.exp(4)])) True """ raise ValueError("A.exp() is only possible if A is a square matrix.")
""" ``A.inv()`` computes the inverse of the square matrix ``A`` (if possible), with the Gauss-Jordan algorithm.
- Raise a ``ValueError`` exception if ``A`` is not square. - Raise a ``ValueError`` exception if ``A`` is singular.
>>> A = Matrix([[1, 2], [3, 4]]) >>> A.inv() [[-2.0, 1.0], [1.5, -0.5]] >>> A * A.inv() == A.inv() * A == eye(A.n) # Rounding mistake can happen (but not here) True >>> Ai = A.inv().round() # No more rounding mistake! >>> A * Ai == Ai * A == eye(A.n) True >>> A.det -2 >>> O = Matrix([[1, 2], [0, 0]]) # O and not 0 >>> O.is_singular True >>> O.inv() # O is singular! Traceback (most recent call last): ... ValueError: A.inv() on a singular matrix (ie. non inversible). >>> O.det 0 """ raise ValueError("A.inv() is only possible if A is a square matrix.") else: except Exception: raise ValueError("A.inv() on a singular matrix (ie. non inversible).")
# ==================================================================== # Gauss elimination process (to get a row echelon form) r""" ``A.gauss()`` implements the Gauss elimination process on matrix ``A``.
When possible, the Gauss elimination process produces a row echelon form by applying linear operations to ``A``.
- If ``maxpivot`` is True, we look for the pivot with higher absolute value (can help reducing rounding mistakes). - If ``verb`` is True, some details are printed at each steps of the algorithm. - ``mode`` can be ``None`` (default), or ``'f'`` for fractions (:class:`Fraction`) or ``'d'`` for decimal numbers (:class:`Decimal`).
- Reference is https://en.wikipedia.org/wiki/Gaussian_elimination#Definitions_and_example_of_algorithm - We chosed to apply rows operations only: it uses elementary operations on lines/rows: :math:`L_i' \longrightarrow L_i - \gamma \times L_k` (method :py:meth:`swap_rows`). - Can swap two columns in order to select the bigger pivot (increases the numerical stability). - The function will raise a ``ValueError`` if the matrix ``A`` is singular (ie. Gauss process cannot conclude). - If ``det`` is ``True``, the returned value is ``c, d`` with ``c`` the row echelon form, and ``d`` the determinant. Reference for this part is `this wikipedia page <https://en.wikipedia.org/wiki/Gaussian_elimination#Computing_determinants>`_.
>>> Matrix([[1, 2], [3, 4]]).gauss() [[1, 2], [0, -2]] >>> Matrix([[1, 2], [1, 2]]).gauss() [[1, 2], [0, 0]] >>> Matrix([[1, 2], [-1, -0.5]]).gauss() [[1, 2], [0, 1.5]] >>> Matrix([[1, 2], [3, 4]]).gauss(maxpivot=True) [[2, 1], [0, 1]] >>> Matrix([[1, 2], [1, 2]]).gauss(maxpivot=True) [[2, 1], [0, 0]] >>> Matrix([[1, 2], [3, 4]]).gauss(det=True) ([[1, 2], [0, -2]], -2) >>> Matrix([[1, 2], [1, 2]]).gauss(det=True) ([[1, 2], [0, 0]], 0) """ # We start with a fresh copy of self.
# Trying to compute the mode ourself mode = 'd' elif mode == 'f': try: currentdet = Fraction(1) c = self.map(Fraction) except Exception as e: print("Failed to convert to Fraction:", e) c = self.copy() elif mode == 'd': try: currentdet = Decimal(1) c = self.map(Decimal) except Exception as e: print("Failed to convert to Decimal:", e) c = self.copy()
i_max = _argmax(list(range(k, m)), [abs(c[k, j]) for j in range(m)]) else: # We found the first i_max such that c[k, i_max] is not zero, or k if none are good # assert c[k, i_max] == max(abs(c[k, j]) for j in range(m))
# XXX Do we really have a singular matrix already ? # raise ValueError("A.gauss_elimination() called on a singular matrix.") # return (c, 0) if det else c # XXX not yet ! # determinant is 0, that is sure at least else: # c.swap_rows(i_max, k) # Swapping two rows multiplies the determinant by -1
print("For the last line, we swapped the {i_max}-th and {k}-th rows, but nothing else.".format(i_max=i_max, k=k))
# Do for all lines below pivot: gamma = Decimal(c[i, k]) / Decimal(c[k, k]) else: # gamma = float(c[i, k]) / float(c[k, k]) # Do for all remaining elements in current line: # We convert to integer if possible, it is prettier :) # FIXME isn't it a cause of rounding mistake? # Fill lower triangular matrix with zeros (because gamma is chosen like that):
# Product of the (-1)**(nb of swaps) * diagonal elements else: # End of gauss()
# ==================================================================== # Gauss-Jordan elimination process (to get a reduced row echelon form) """ ``A.gauss_jordan()`` implements the Gauss elimination process on matrix ``A``.
- If ``inv`` is ``True``, the returned value is ``J_n, A**(-1)`` with ``J_n`` the reduced row echelon form of ``A``, and ``A**(-1)`` the computed inverse of A. - If ``maxpivot`` is ``True``, we look for the pivot with higher absolute value (can help reducing rounding mistakes). """ # We start with a fresh copy of self.
# Trying to compute the mode ourself elif all(isinstance(x, Decimal) for x in self): mode = 'd' except Exception as e: print("Failed to convert to Fraction:", e) c = self.copy() elif mode == 'd': try: c = self.map(Decimal) except Exception as e: print("Failed to convert to Decimal:", e) c = self.copy()
raise ValueError("A.gauss_jordan(inv=True) is only possible if A is a square matrix.") # print("OK, trying to produce the inverse of c also.") # print("At first, c is:", c) # print("Experimental construction of the augmented matrix of size (n, 2m) (with n = {}, m = {}).".format(c.n, c.m)) # c = mat_from_f(lambda i, j: c[i, j] if j < c.m else int((i+c.m)==j), c.n, 2 * c.m) # # Reference is https://en.wikipedia.org/wiki/Augmented_matrix # print("Now, c is:", c) print("\nj =", j, "=> current c is:", c) if inv: print("INFO: Current cinv is:", cinv) print("Looking for the pivot with r = {r}, j = {j}.".format(r=r + 1, j=j)) print("Indeces:", list(range(r + 1, n))) print("Values:", [abs(c[i, j]) for i in range(n)]) # k = _argmax(range(r+1, n), [abs(c[i, j]) for i in range(n)]) k = _argmax(list(range(r + 1, n)), [abs(c[i, j]) for i in range(n)]) else: for _ in range(r + 1, n): if c[k, j] != 0: k = k break # We found the first i_max such that c[k, i_max] is not zero, or k if none are good print("For r = {}, c[k, j] is the pivot (k = {}, j = {}), equals to {}.".format(r + 1, k, j, c[k, j])) raise ValueError("A.gauss_jordan() called on a singular matrix.") print("Pivot is not zero ({}), so we divide the k-th (k = {}) row by the pivot.".format(pivot, k)) # c[k, :] /= pivot # Divising the k-th row by the pivot # c[k, jjjj] /= pivot elif mode == 'd': c[k, jjjj] = Decimal(c[k, jjjj] / pivot) else: c[k, jjjj] /= pivot print("INFO: Same linear operation is applied to cinv: cinv[k, :] /= pivot") # cinv[k, :] /= pivot # cinv[k, jjjj] /= pivot elif mode == 'd': cinv[k, jjjj] = Decimal(cinv[k, jjjj] / pivot) else: cinv[k, jjjj] /= pivot if verb: print("We swap the rows r = {} and k = {}.".format(r, k)) print("Before: R_k =", c[k, :], "R_r =", c[r, :]) c.swap_rows(k, r) # Swap the rows r and k # c.swap_cols(k, r) # Swap the columns r and k if inv: if verb: print("INFO: Same linear operation is applied to cinv: cinv.swap_rows(k, r)") cinv.swap_rows(k, r) # Swap the rows r and k # cinv.swap_cols(k, r) # Swap the columns r and k if verb: print("After: R_k =", c[k, :], "R_r =", c[r, :]) print("For i = {}.".format(i)) print("Before: R_{i} =", c[i, :]) print("R_{i} <-- R_{i} - c[{i}, {j}] * R_{r}. c[{i}, {j}] is {cij}".format(i=i, j=j, r=r, cij=cij)) print("INFO: Same linear operation is applied to cinv: cinv[i, :] -= cinv[r, :] * c[i, j]") # R_i <-- R_i - c[i,j] * R_r print("After: R_{i} =", c[i, :]) else: print("No modification here (i = r or c[i, j] = 0).") # Done print("Done ! c =", c) # c = c.map(lambda x: int(x) if int(x) == x else x) # Pretty priting only print("Done ! cinv =", cinv) else: return c
# ==================================================================== # Applications of the Gauss elimination process
""" ``A.rank`` uses the Gauss elimination process to compute the rank of the matrix ``A``, by simply counting the number of non-zero elements on the diagonal of the echelon form.
.. todo:: The Gauss process (:py:meth:`gauss`) has to be changed, and improved for singular matrices (when the rank is not maximum!).
>>> Matrix([[1, 2], [3, 4]]).rank 2 >>> Matrix([[1, 2], [1, 2]]).rank 1 >>> zeros(7).rank 0 >>> eye(19).rank 19 """
""" ``A.det`` uses the Gauss elimination process to compute the determinant of the matrix ``A``.
.. note:: Because it depends of the number of elementary operations performed in the Gauss method, we had to modify the :py:meth:`gauss` method...
>>> Matrix([[1, 2], [3, 4]]).det -2 >>> Matrix([[1, 2], [1, 2]]).det 0 >>> zeros(7).det 0 >>> eye(19).det 1 """
# ==================================================================== # Extra methods (not required in the project)
""" ``A.count(value)`` counts how many times the element ``value`` is in the matrix ``A``.
>>> Matrix([[1, 2], [3, 4]]).count(2) 1 >>> Matrix([[1, 2], [1, 2]]).count(2) 2 >>> zeros(7).count(2) 0 >>> zeros(7).count(0) 49 >>> eye(19).count(1) 19 >>> eye(19).count(0) 342 """
""" ``value in A`` <-> ``A.__contains__(value)`` tells if the element ``value`` is present in the matrix ``A``.
>>> 4 in Matrix([[1, 2], [3, 4]]) True >>> 4 in Matrix([[1, 2], [1, 2]]) False >>> O, I = zeros(7), eye(7) >>> 3 * I**2 + 2 * I + O ** 0 [[6, 0, 0, 0, 0, 0, 0], [0, 6, 0, 0, 0, 0, 0], [0, 0, 6, 0, 0, 0, 0], [0, 0, 0, 6, 0, 0, 0], [0, 0, 0, 0, 6, 0, 0], [0, 0, 0, 0, 0, 6, 0], [0, 0, 0, 0, 0, 0, 6]] >>> 6 in (3 * I**2 + 2 * I + O ** 0) True """
""" Apply the function ``f`` to each of the coefficient of the matrix ``A`` (returns a new matrix).
>>> O, I = zeros(2), eye(2) >>> I.map(lambda x: x * 4) [[4, 0], [0, 4]] >>> O.map(lambda x: x + 6) [[6, 6], [6, 6]] >>> A = Matrix([[-1j, -2j], [-2j, -1j]]) >>> A.map(lambda z: abs(z)) [[1.0, 2.0], [2.0, 1.0]] >>> A.map(lambda z: int(abs(z))) [[1, 2], [2, 1]] >>> A.map(lambda z: z + 1j) [[0j, -1j], [-1j, 0j]] >>> A.map(lambda z: '"%s"' % str(z)) [["-1j", "-2j"], ["-2j", "-1j"]] >>> A.map(lambda z: "Look: %s" % str(z)) [[Look: -1j, Look: -2j], [Look: -2j, Look: -1j]]
- If ``f`` needs arguments or key-words arguments, use the ``*args`` and ``**kwargs`` :
>>> def f(x, n, offset=0): ... return (x ** n) + offset >>> A = Matrix([[1, 2], [2, 1]]) >>> A.map(f, 2) [[1, 4], [4, 1]] >>> A.map(f, 2, offset=4) [[5, 8], [8, 5]] """
""" ``A.round([ndigits=8])`` <-> rounds every coefficient of ``A`` to ``ndigits`` digits after the comma.
>>> A = (1. / 3.) * eye(2) + 4 >>> A.round(0) [[4.0, 4.0], [4.0, 4.0]] >>> A.round(2) [[4.33, 4.0], [4.0, 4.33]] >>> A.round(7) [[4.3333333, 4.0], [4.0, 4.3333333]] """ return self.map(lambda x: round(x, ndigits))
# ==================================================================== # Iterating over a matrix
""" ``iter(A)`` <-> ``A.__iter__()`` is used to create an iterator from the matrix ``A``.
- The values are looped rows by rows, then columns then columns. - This method is called when an iterator is required for a container. This method should return a new iterator object that can iterate over all the objects in the container.
>>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> list(A) [1, 2, 3, 4, 5, 6, 7, 8, 9] """
# next and __next__ for iterating over the values of our matrix
""" For Python 3 compatibility.""" return self.next()
""" Generator for iterating the matrix ``A``.
- The values are looped rows by rows, then columns then columns.
>>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> for x in A: ... print(x) 1 2 3 4 5 6 7 8 9 >>> for i, x in enumerate(A): ... print(i, "th value of A is", x) 0 th value of A is 1 1 th value of A is 2 2 th value of A is 3 3 th value of A is 4 4 th value of A is 5 5 th value of A is 6 6 th value of A is 7 7 th value of A is 8 8 th value of A is 9 """ if (self._i0 < self.n) and (self._j0 < self.m): v = self[self._i0, self._j0] if self._i0 < self.n - 1: self._i0 += 1 else: self._i0 = 0 self._j0 += 1 return v else: raise StopIteration()
# ==================================================================== # To deal nicely with matrices of complex numbers
def real(self): """ Real part of the matrix ``A``, coefficient wise.
>>> A = Matrix([[1j, 2j], [2j, 1j]]) >>> A.real [[0.0, 0.0], [0.0, 0.0]] >>> A = Matrix([[1+6j, 2], [-1+2j, 1+9j]]) >>> A.real [[1.0, 2], [-1.0, 1.0]] """
def imag(self): """ Imaginary part of the matrix ``A``, coefficient wise.
>>> A = Matrix([[-1j, -2j], [-2j, -1j]]) >>> A.imag [[-1.0, -2.0], [-2.0, -1.0]] """
""" Conjugate part of the matrix ``A``, coefficient wise.
>>> A = Matrix([[-1j, -2j], [-2j, -1j]]) >>> A.conjugate() [[1j, 2j], [2j, 1j]] """
# ==================================================================== # Dot product and norm
r""" ``A.dot(v)`` computes the dot multiplication of the matrix ``A`` and the vector ``v`` (:math:`A \dot v`).
- ``v`` can be a matrix (:class:`Matrix`) of size ``(m, 1)``, or a list of size ``m``.
>>> A = Matrix([[1, 1], [1, -1]]) >>> v = [2, 3] >>> A.dot(v) [[5], [-1]] >>> v = Matrix([[2], [-3]]) >>> A.dot(v) [[-1], [5]]
.. warning:: An exception ``ValueError`` is raised if the sizes does not allow the dot product:
>>> A.dot(v.T) # v.T is not a column vector! Traceback (most recent call last): ... ValueError: A.dot(v): the vector v = [[2, -3]] is not a vector: v.m = 2 != 1. >>> v = Matrix([[2], [-3], [7]]) >>> A.dot(v) Traceback (most recent call last): ... ValueError: A.dot(v): the size of the vector v = [[2], [-3], [7]] should be compatible with the size of the matrix self = [[1, 1], [1, -1]]. Here self.m = 2 and v.n = 3, are different. >>> v = [1, 2, 3, 4, 5] >>> A.dot(v) Traceback (most recent call last): ... ValueError: A.dot(v): the size of the vector v = [[1], [2], [3], [4], [5]] should be compatible with the size of the matrix self = [[1, 1], [1, -1]]. Here self.m = 2 and v.n = 5, are different. """ elif v.m != 1: raise ValueError(("A.dot(v): the vector v = {} is not a vector: v.m = {} != 1.".format(v, v.m))) elif self.m != v.n: raise ValueError(("A.dot(v): the size of the vector v = {} should be compatible with the size of the matrix self = {}. Here self.m = {} and v.n = {}, are different.".format(v, self, self.m, v.n))) else: # Convert the iterator v to a list, then to a column vector except Exception: raise ValueError(("A.dot(v): impossible to convert v = {} to a column vector.".format(v)))
r""" ``A.norm(p)`` computes the p-norm of the matrix ``A``, default is ``p = 2``.
- Mathematically defined as p-root of the sum of the p-power of *modulus* of its coefficients :
.. math:: \|A\|_{p} := \left( \sum\limits_{1 \leq i \leq n, 1 \leq j \leq m} {|A_{i,j}|}^p \right)^{\frac{1}{p}}
- If ``p = 'inf'``, the max norm is returned (ie. infinity norm), defined by :math:`\|A\|_{\infty} := \max_{i,j} |A_{i,j}|`. - Reference is `Matrix norm (on Wikipedia) <https://en.wikipedia.org/wiki/Matrix_norm#.22Entrywise.22_norms>`_.
>>> A = Matrix([[1, 2], [-3, -1]]) >>> A.norm() # (1)**2 + (2)**2 + (-3)**2 + (-1)**2 3.872983346207417 >>> 15**0.5 3.872983346207417 >>> A.norm('inf') 3 >>> A.norm(1) == 7 # (1) + (2) + (3) + (1) True >>> A.norm(3) 3.332221851645953 """ return max(abs(x) for x in self) else:
""" ``A.normalized()`` return a new matrix, which **columns vectors are normalized** by using the norm ``2`` (or the given function ``fnorm``).
- Will **not fail** if a vector has norm ``0`` (it is just not modified). - Reference is `Orthogonalization (on Wikipedia) <https://en.wikipedia.org/wiki/Orthogonalization>`_. - Any extra arguments ``args``, ``kwargs`` are given to the function ``fnorm``.
>>> A = Matrix([[1, 2], [-3, -1]]) >>> A.normalized(p='inf') # doctest: +ELLIPSIS [[0.333333..., 1.0], [-1.0, -0.5]] >>> eye(5).normalized(p='inf').map(int) # normalize then round to an int [[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1]] >>> B = -eye(5) >>> (2*B).normalized() # each vector is divided by its norm = 2 [[-1.0, 0.0, 0.0, 0.0, 0.0], [0.0, -1.0, 0.0, 0.0, 0.0], [0.0, 0.0, -1.0, 0.0, 0.0], [0.0, 0.0, 0.0, -1.0, 0.0], [0.0, 0.0, 0.0, 0.0, -1.0]] >>> B.normalized(p='inf') [[-1.0, 0.0, 0.0, 0.0, 0.0], [0.0, -1.0, 0.0, 0.0, 0.0], [0.0, 0.0, -1.0, 0.0, 0.0], [0.0, 0.0, 0.0, -1.0, 0.0], [0.0, 0.0, 0.0, 0.0, -1.0]]
It works also for a simple vector:
>>> v = Matrix([[1], [-2], [3]]) >>> v.normalized() # doctest: +ELLIPSIS [[0.267261...], [-0.534522...], [0.801783...]] >>> v.normalized(p=2) # doctest: +ELLIPSIS [[0.267261...], [-0.534522...], [0.801783...]] >>> v.normalized() * (14**0.5) [[1.0], [-2.0], [3.0]] >>> v.normalized(p=1) # doctest: +ELLIPSIS [[0.166666...], [-0.333333...], [0.5]] >>> v.normalized(p=1) * 6 [[1.0], [-2.0], [3.0]] >>> 6 * v.normalized(p=1) [[1.0], [-2.0], [3.0]] """ if fnorm is None: def fnorm(x, *args, **kwargs): return x.norm(*args, **kwargs) c = self.copy() for j in range(c.cols): normofthiscol = fnorm(c[:, j], *args, **kwargs) if normofthiscol != 0: for i in range(c.rows): c[i, j] /= normofthiscol return c
""" ``abs(A)`` <-> ``A.abs()`` <-> ``A.__abs__()`` computes the absolute value / modulus of ``A`` coefficient-wise.
>>> A = Matrix([[-4, 2+2j], [0, 4j]]) >>> abs(A) # doctest: +ELLIPSIS [[4, 2.828427...], [0, 4.0]] >>> B = -eye(2) >>> B.abs() [[1, 0], [0, 1]] """
# ==================================================================== # Trace and other special values
r""" ``A.trace()`` computes the trace of ``A`` :
.. math:: \mathrm{Tr}(A) := \sum\limits_{1 \leq i \leq \min(n, m)} A_{i, i}
>>> A = Matrix([[-4, 2+2j], [0, 4j]]) >>> A.trace() (-4+4j) >>> eye(19).trace() 19 >>> zeros(20).trace() 0 >>> ones(100).trace() 100 """
# TODO: Try to find an algorithm to approximatively compute eigen values, and eigen vectors ?
# ==================================================================== # Check if a matrix is square, symetric, anti-symetric, diagonal, skew-symetric (hermitian) etc
def is_square(self): """ ``A.is_square`` tests if ``A`` is **square** or not.
>>> A = Matrix([[-4, 2+2j], [0, 4j]]) >>> A.is_square True >>> v = Matrix([[-4], [0]]) >>> v.is_square False """
def is_symetric(self): """ ``A.is_symetric`` tests if ``A`` is **symetric** or not.
>>> A = Matrix([[-4, 2+2j], [0, 4j]]) >>> A.is_symetric False >>> eye(30).is_symetric True """
def is_anti_symetric(self): """ ``A.is_anti_symetric`` tests if ``A`` is **anti-symetric** or not.
>>> A = Matrix([[0, 1], [-1, 0]]) >>> A.is_anti_symetric True >>> eye(30).is_anti_symetric False """
def is_diagonal(self): """ ``A.is_diagonal`` tests if A is **diagonal** or not.
>>> eye(40).is_diagonal True >>> A = Matrix([[0, 1], [-1, 0]]) >>> A.is_diagonal False >>> A = diag(range(30)) >>> A.is_diagonal True """
def is_hermitian(self): r""" ``A.is_hermitian`` tests if ``A`` is **Hermitian** or not (tests if :math:`A^{*} = A`, ie. ``conjugate(A.T) == A)``).
>>> A = Matrix([[1, 2j], [-2j, 1]]) >>> A.is_hermitian True >>> eye(30).is_hermitian True >>> (1j * ones(3)).is_hermitian False """
def is_lower(self): """ ``A.is_lower`` tests if ``A`` is **lower triangular** or not.
>>> A = Matrix([[8, 1], [0, 7]]) >>> A.is_lower False >>> A.T.is_lower True """
def is_upper(self): """ ``A.is_upper`` tests if ``A`` is **upper triangular** or not.
>>> A = Matrix([[2, 0], [3, 4]]) >>> A.is_upper False >>> A.T.is_upper True """
def is_zero(self): """ ``A.is_zero`` tests if ``A`` is the **zero matrix** or not.
>>> A = Matrix([[2, 0], [3, 4]]) >>> A.is_zero False >>> zeros(30).is_zero True >>> (0 * A).is_zero True """
def is_singular(self): """ ``A.is_singular`` tests if ``A`` is **singular** (ie. non-invertible) or not.
.. note:: It computes the determinant by using the Gauss elimination process (:py:meth:`det`).
>>> A = Matrix([[2, 0], [3, 4]]) >>> A.is_singular False >>> zeros(3).is_singular True >>> (0 * A).is_singular True >>> Matrix([[2, 0], [4, 0]]).is_singular True """ return self.det == 0
# ==================================================================== # Linear operations *in place*
""" ``A.swap_cols(j1, j2)`` changes *in place* the ``j1``-th and ``j2``-th *columns* of the matrix ``A``.
>>> A = Matrix([[2, 0], [3, 4]]); A [[2, 0], [3, 4]] >>> A.swap_cols(0, 1); A [[0, 2], [4, 3]] """
""" ``A.swap_rows(i1, i2)`` changes *in place* the ``i1``-th and ``i2``-th *rows* of the matrix ``A``.
>>> A = Matrix([[2, 0], [3, 4]]); A [[2, 0], [3, 4]] >>> A.swap_rows(0, 1); A [[3, 4], [2, 0]] """
# ======================================================================== # Adjugate matrix (https://en.wikipedia.org/wiki/Adjugate_matrix)
r""" ``A.minor(i, j)`` <-> ``minor(A, i, j)`` returns the ``(i, j)`` minor of ``A``, defined as the determinant of the submatrix ``A[i0, j0]`` for ``i0 != i`` and ``j0 != j``.
- Complexities: memory is :math:`\mathcal{O}(n^2)`, time is :math:`\mathcal{O}(n^3)` (1 determinant of size ``n - 1``).
>>> A = Matrix([[1, 2], [3, 4]]) >>> A.minor(0, 0) 4 >>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> A.minor(0, 0) # | 5 6 8 9 | = 5 * 9 - 6 * 8 = -3 -3.000000000000007 >>> A.minor(1, 0) # | 2 3 8 9 | = 2 * 9 - 3 * 8 = -6 -6 """
r""" ``A.cofactor(i, j)`` <-> ``cofactor(A, i, j)`` returns the ``(i, j)`` cofactor of ``A``, defined as the ``(-1)**(i + j)`` times to ``(i, j)`` minor of ``A`` (cf. :py:meth:`minor`).
- Complexities: memory is :math:`\mathcal{O}(n^2)`, time is :math:`\mathcal{O}(n^3)` (1 determinant of size ``n - 1``).
>>> A = Matrix([[1, 2], [3, 4]]) >>> A.cofactor(0, 0) 4 >>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> A.cofactor(0, 0) # (-1)**0 * | 5 6 8 9 | = 5 * 9 - 6 * 8 = -3 -3.000000000000007 >>> A.cofactor(1, 0) # (-1)**1 * | 2 3 8 9 | = -(2 * 9 - 3 * 8) = 6 6 """
r""" ``A.adjugate()`` <-> ``adjugate(A)`` returns the **adjugate matrix** of ``A``.
- Reference is https://en.wikipedia.org/wiki/Adjugate_matrix#Inverses. - Complexities: memory is :math:`\mathcal{O}(n^2)`, time is :math:`\mathcal{O}(n^5)` (:math:`n^2` determinants of size ``n - 1``). - Using the adjugate matrix for computing the inverse is a BAD method : too time-consuming ! LU or Gauss-elimination is only :math:`\mathcal{O}(n^3)`.
>>> A = Matrix([[2, 0], [3, 4]]) >>> A.adjugate() [[4, -3], [0, 2]] >>> A * A.adjugate() == A.det * eye(A.n) False >>> A * A.adjugate().T == A.det * eye(A.n) True """
""" ``A.type()`` returns the matrix of types of coefficients of ``A``. """
# End of that class Matrix # ========================================================================
# ======================================================================== # Utility functions
""" ``ones(n, m)`` is a matrix of size ``(n, m)`` filled with ``1``.
>>> ones(3, 2) [[1, 1], [1, 1], [1, 1]] >>> ones(2, 3) [[1, 1, 1], [1, 1, 1]]
- It works with only one dimension, or with a tuple ``(n, m)`` :
>>> ones(2) [[1, 1], [1, 1]] >>> ones((2, 3)) [[1, 1, 1], [1, 1, 1]] """ n, m = n
""" ``zeros(n, m)`` is a matrix of size ``(n, m)`` filled with ``0``.
>>> zeros(3, 2) [[0, 0], [0, 0], [0, 0]] >>> zeros(2, 3) [[0, 0, 0], [0, 0, 0]] >>> ones(2, 3) == zeros(2, 3) + 1 True >>> zeros(2, 3) == ones(2, 3) * 0 True
- It works with only one dimension, or with a tuple ``(n, m)`` :
>>> zeros(2) [[0, 0], [0, 0]] >>> zeros((2, 3)) [[0, 0, 0], [0, 0, 0]] """ n, m = n
""" ``eye(n)`` is the (square) identity matrix of size ``(n, n)`` (``1`` on the diagonal, ``0`` outside).
>>> eye(2) [[1, 0], [0, 1]] >>> zeros(18) == eye(18) * 0 True >>> eye(60).is_diagonal True >>> eye(40).is_square True >>> eye(20).is_singular False >>> eye(5).det 1 >>> eye(7).trace() 7 """
""" ``diag(d)`` creates a matrix from a list ``d`` (or iterator) of diagonal values, or with ``n``-times the value ``d`` if ``d`` is not an iterator and ``n`` is an integer.
>>> D = diag(range(1,6)) >>> D[2, :] [[0, 0, 3, 0, 0]]
We can check the usual properties of diagonal matrices:
>>> D.trace() 15 >>> D.trace() == sum(range(1,6)) True >>> D.det 120 >>> from math import factorial >>> D.det == factorial(5) True
Other examples:
>>> diag([-1, 1]) [[-1, 0], [0, 1]] >>> diag([-4, 1]) + 3 [[-1, 3], [3, 4]]
We can also use the optional argument ``n``:
>>> diag(3.14, 3) [[3.14, 0, 0], [0, 3.14, 0], [0, 0, 3.14]] >>> diag([3.14]*3) # Same ! [[3.14, 0, 0], [0, 3.14, 0], [0, 0, 3.14]] """ return diag([d] * n) else:
""" ``mat_from_f(f, n, m=None)`` creates a matrix of size ``(n, m)`` initialized with the function ``f`` : ``A[i, j] = f(i, j)``.
- Default value for ``m`` is ``n`` (square matrix).
.. warning:: ``f`` has to accept (at least) two arguments ``i, j``.
>>> mat_from_f(lambda i, j: 1 if i == j else 0, 3) == eye(3) True >>> mat_from_f(lambda i, j: 1, 3) == ones(3) True >>> mat_from_f(lambda i, j: i+j, 3) [[0, 1, 2], [1, 2, 3], [2, 3, 4]] >>> mat_from_f(lambda i, j: i*j, 3) [[0, 0, 0], [0, 1, 2], [0, 2, 4]]
- Any extra arguments ``args``, ``kwargs`` are given to the function ``f``.
>>> def f(i, j, e, offset=0): ... return (i * e) + offset >>> mat_from_f(f, 2, 2, 4) # n = 2, m = 2, e = 4 [[0, 0], [4, 4]] >>> mat_from_f(f, 2, 2, 4, offset=10) # n = 2, m = 2, e = 4, offset = 10 [[10, 10], [14, 14]]
- Remark: it is similar to ``Array.make`` (or ``Array.init``) in OCaml (v3.12+) or ``String.create`` (or ``String.make``). """
# ======================================================================== # Functions that are just wrappers around methods
r""" ``det(A)`` <-> ``A.det`` computes the determinant of ``A`` (in :math:`\mathcal{O}(n^3)`).
>>> det(eye(2)) 1 >>> det((-1) * eye(4)) 1 >>> det((-1) * eye(5)) -1 """
r""" ``rank(A)`` <-> ``A.rank`` computes the rank of ``A`` (in :math:`\mathcal{O}(n^3)`).
>>> rank(eye(2)) 2 """
r""" ``gauss(A)`` <-> ``A.gauss()`` applies the Gauss elimination process on ``A`` (in :math:`\mathcal{O}(n^3)`). """ return A.gauss_elimination(*args, **kwargs)
r""" ``gauss_jordan(A)`` <-> ``A.gauss_jordan()`` applies the Gauss-Jordan elimination process on ``A`` (in :math:`\mathcal{O}(n^3)`). """ return A.gauss_jordan(*args, **kwargs)
r""" ``inv(A)`` <-> ``A.inv()`` **tries** to compute the inverse of ``A`` (in :math:`\mathcal{O}(n^3)`).
>>> inv(eye(2)) == eye(2) True """
r""" ``exp(A)`` <-> ``A.exp()`` computes an approximation of the exponential of ``A`` (in :math:`\mathcal{O}(n^3 * limit)`).
>>> import math >>> e = math.exp(1.0) >>> C = diag([1, 4]) >>> exp(C) == diag([e ** 1, e ** 4]) == diag([math.exp(1), math.exp(4)]) # Rounding mistakes! False >>> exp(C).almosteq(diag([e ** 1, e ** 4])) # No more rounding mistakes! True >>> diag([e ** 1, e ** 4]).almosteq(diag([math.exp(1), math.exp(4)])) True """
# ======================================================================== # The LU decomposition
r""" ``PLUdecomposition(A)`` computes the **permuted LU decomposition** for the matrix ``A``.
- Operates in time complexity of :math:`\mathcal{O}(n^3)`, memory of :math:`\mathcal{O}(n^2)`. - ``mode`` can be ``None`` (default), or ``'f'`` for fractions (:class:`Fractions`) or ``'d'`` for decimal (:class:`Decimal`) numbers. - Returned ``P, L, U`` that satisfies ``P*A = L*U``, with ``P`` being a permutation matrix, ``L`` a lower triangular matrix, ``U`` an upper triangular matrix. - Will raise a ``ValueError`` exception if ``A`` is singular. - Reference is `Gauss elimination (on Wikipedia) <https://en.wikipedia.org/wiki/Gaussian_elimination#Definitions_and_example_of_algorithm>`_. - We chosed to apply rows operations only: it uses elementary operations on lines/rows: :math:`L_i' \longrightarrow L_i - \gamma \times L_k` (method swap_rows). - Can swap two columns in order to select the bigger pivot (increases the numerical stability). """
# Trying to compute the mode ourself elif all(isinstance(x, Decimal) for x in U): mode = 'd' elif mode == 'f': try: U = U.map(Fraction) L = L.map(Fraction) except Exception as e: print("Failed to convert to Fraction:", e) U = U.copy() elif mode == 'd': try: U = U.map(Decimal) L = L.map(Decimal) except Exception as e: print("Failed to convert to Decimal:", e) U = U.copy()
# Now we can start # Find the pivot on the k-th row # Is the pivot U[k, k] is zero, we find a possible better pivot # Is the pivot U[k, i_max] is still non zero, we cannot do anything, because the matrix is singular. raise ValueError("PLUdecomposition() has been called on a singular matrix.") else: # U.swap_cols(i_max, k) # The matrix P will keep track of the permutations performed during the Gaussian Elimination process: # P.swap_cols(i_max, k)
# Do for all lines/rows below pivot: elif mode == 'd': gamma = Decimal(U[i, k]) / Decimal(U[k, k]) else: # gamma = float(U[i, k]) / float(U[k, k]) gamma = U[i, k] / U[k, k] # Do for all remaining elements in current line: # Add - gamma times row k to row i of U # We convert to integer if possible, it is prettier :) # if int(U[i, j]) == U[i, j]: # U[i, j] = int(U[i, j]) # Fill lower triangular matrix with zeros (because gamma is chosen like that): elif mode == 'd': U[i, k] = Decimal(0) else: U[i, k] = 0 # The entries of L below the diagonal are gradually replaced by the negatives of multiples used in the corresponding row operations of type #1.
# Moreover, any pair of entries that both lie below the diagonal # in these same two rows (i_max and k) of L must also be interchanged, # while entries lying on and above its diagonal need to stay in their place. # L[j, i_max], L[j, k] = L[j, k], L[j, i_max]
# End of PLUdecomposition()
# ======================================================================== # Other functions
""" ``norm(A, p)`` <-> ``A.norm(p)`` computes the p-norm of ``A`` (default is ``p = 2``)."""
""" ``trace(A)`` <-> ``A.trace()`` computes the trace of ``A``."""
# ======================================================================== # Random matrix
""" ``rand_matrix(n, m, k)`` generates a new random matrix of size ``(n, m)`` with each coefficients being integers, randomly taken between ``-k`` and ``k`` (bound *included*).
>>> from random import seed >>> seed(0) # We want the examples to always be the same >>> rand_matrix(2, 3) [[7, 5, -2], [-5, 0, -2]] >>> rand_matrix(3, 2, 40) [[23, -16], [-2, 7], [33, 0]] >>> rand_matrix(4, 4, 100) [[-44, 51, 24, -50], [82, 97, 62, 81], [-38, 46, 80, 37], [-6, -80, -13, 22]] """
""" ``rand_matrix_float(n, m, k)`` generates a new random matrix of size ``(n, m)`` with each coefficients being float numbers, randomly taken between ``-k`` and ``k`` (right bound excluded).
>>> from random import seed >>> seed(0) # We want the examples to always be the same >>> rand_matrix_float(2, 3) [[6.8884370305, 5.15908805881, -1.58856838338], [-4.82166499414, 0.225494427372, -1.90131725099]] >>> rand_matrix_float(3, 2, 1) [[0.56759717807, -0.393374547842], [-0.0468060916953, 0.16676407891], [0.816225770391, 0.00937371163478]] >>> rand_matrix_float(4, 4, 20) [[-8.72648622401, 10.2321681663, 4.73475986701, -9.9797463455], [16.3898502387, 19.3114190415, 12.4086894399, 16.0866380176], [-7.59409722723, 9.19326993041, 15.9535315187, 7.35935727662], [-1.11429138189, -15.9719516773, -2.63312658185, 4.43547893775]] """
# ======================================================================== # 2 auxiliary functions used by the Gauss elimination process
""" Compute the index ``i`` in ``indexes`` such that the ``array[i]`` is the bigger.""" raise ValueError("argmax() arg indexes is a non-valid sequence.") # bestvalue = array[indexes[0]] raise ValueError("argmax() arg is a non-valid sequence.")
""" Compute the product of the values in the iterator ``iterator``. Empty product is 1."""
# 2 auxiliary functions for implementing the generalized __setitem__ method
""" ``b if (a is None), else a``.
- Useful for converting a ``slice`` object to a ``range`` object (:class:`slice`, :class:`range`). """ return b if (a is None) else a
""" Get a ``range`` of indeces from a ``slice`` object (:class:`slice`, :class:`range`).
- Thanks to `this answer on stackoverflow.com <http://stackoverflow.com/a/13855369>`_. """ return range(_ifnone(sliceobject.start, 0), sliceobject.stop, _ifnone(sliceobject.step, 1))
# ======================================================================== # The Gram-Schmidt orthogonalization process
r""" (Hermitian) dot product of the two vectors ``vx`` and ``vy`` (sum of ``conjugate(vx[i]) * vy[i]``) :
.. math:: \mathbf{x} . \mathbf{y} = \langle \mathbf{x}, \mathbf{y} \rangle := \sum_{1 \leq i \leq n} \overline{x_i} \times y_i.
>>> vx = [1, 2, 3]; vy = [-1, 0, 4] >>> innerproduct(vx, vy) 11
.. warning:: The conjugate is on the first vector, as always for Hermite spaces and Hermitian inner product.
>>> vx = [1j, 2j, 3j]; vy = [-1, 0, 4] >>> (-1j) * (-1) + (-2j) * (0) + (-3j) * (4) -11j >>> innerproduct(vx, vy) -11j """ assert len(vx) == len(vy) res = 0 # XXX Typo in the subject for x, y in zip(vx, vy): if hasattr(x, "conjugate"): res += x.conjugate() * y else: res += x * y return res # sum(x.conjugate() * y for x, y in zip(vx, vy))
r""" Shortcut for the square of the norm of the vector ``u``:
.. math:: \| u \|^2 := \langle u, u \rangle.
>>> u = [1, 2, 3] >>> norm_square(u) 14
- It works for imaginary valued vectors:
>>> u = [1j, -2j, 3j] >>> norm_square(u) 14.0
- And it also works for complex valued vectors:
>>> u = [1+1j, 2-2j, 3+3j] >>> norm_square(u) 28.0 """ res = innerproduct(u, u) if hasattr(res, "real"): return res.real else: return res
r""" Shortcut for the canonical norm of the vector ``u``:
.. math: \| u \| = \sqrt{\langle u, u \rangle}.
>>> u = [1, 2, 3] >>> norm2(u) 3.7416573867739413
- It works for imaginary valued vectors:
>>> u = [1j, -2j, 3j] >>> norm2(u) 3.7416573867739413
- And it also works for complex valued vectors:
>>> u = [1+1j, 2-2j, 3+3j] >>> norm2(u) 5.291502622129181 """ return norm_square(u) ** 0.5
""" Multiply the vector ``vx`` by the constant ``c`` (scalar, ie. real or complex).
>>> vx = [1, 2, 3]; vy = [-1, 0, 4] >>> vect_const_multi(vx, 2) [2, 4, 6] >>> vect_const_multi(vy, -4) [4, 0, -16] """ return [c * x for x in vx]
r""" Projection of the vector ``v`` into the vector ``u`` (:math:`\mathrm{proj}_u(v)` as written on Wikipedia).
>>> u = [1, 2, 3]; v = [-1, 0, 4] >>> proj(u, v) # 11/14 * u [0.7857142857142857, 1.5714285714285714, 2.357142857142857] >>> proj(u, v) == [(11/14) * x for x in u] True """ nsqu = norm_square(u) if nsqu == 0: return [0] * len(u) else: # udotu = float(nsqu) # useless, I imported division from __future__ return vect_const_multi(u, innerproduct(u, v) / nsqu)
r""" Basic implementation of the Gram-Schmidt process for the column vectors of the matrix ``V``, in the easy case of :math:`\mathbb{R}^n` with the usual dot product.
- The matrix is interpreted as a family of *column* vectors. - Reference for notations, concept and proof is `Gram-Schmidt process (on Wikipedia) <https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process>`_. - If ``normalized`` is ``True``, the vectors are normalized before being returned.
>>> V = Matrix([[1, 2, 3], [-1, 0, 4]]) >>> gram_schmidt(V) [[1, 2, 3], [-1, 0, 4]] """ n, m = V.n, V.m U = V.copy()
for k in range(1, n): # U[k, :] -= sum_j(proj(U[j, :], U[k, :])) for j in range(0, k - 1): p = proj(U[j, :], U[k, :]) for t in range(m): U[k, t] -= p[t] # Now u_{k} is orthogonal with all the previous u_{j} (j < k) ! if normalized: return U.normalized() else: return U
# ======================================================================== # Adjugate matrix (https://en.wikipedia.org/wiki/Adjugate_matrix)
r""" ``minor(A, i, j)`` <-> ``A.minor(i, j)`` returns the ``(i, j)`` minor of ``A``, defined as the determinant of the submatrix ``A[i0,j0]`` for ``i0 != i`` and ``j0 != j``.
- Complexities: memory is :math:`\mathcal{O}(n^2)`, time is :math:`\mathcal{O}(n^3)` (1 determinant of size ``n - 1``). """
r""" ``cofactor(A, i, j)`` <-> ``A.cofactor(i, j)`` returns the ``(i, j)`` cofactor of ``A``, defined as the ``(-1)**(i + j)`` times to ``(i, j)`` minor of ``A`` (cf. :py:func:`minor`).
- Complexities: memory is :math:`\mathcal{O}(n^2)`, time is :math:`\mathcal{O}(n^3)` (1 determinant of size ``n - 1``). """
r""" ``adjugate(A)`` <-> ``A.adjugate()`` returns the adjugate matrix of ``A``.
- Reference is `Adjugate matrix (on Wikipedia) <https://en.wikipedia.org/wiki/Adjugate_matrix#Inverses>`_. - Complexities: memory is :math:`\mathcal{O}(n^2)`, time is :math:`\mathcal{O}(n^5)` (n^2 determinants of size ``n - 1``). - Using the adjugate matrix for computing the inverse is a BAD method : too time-consuming ! LU or Gauss-elimination is only :math:`\mathcal{O}(n^3)`. """
# ======================================================================== # TODO Solver for linear equation A.x = b # Use the pseudo-inverse ? https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse#The_iterative_method_of_Ben-Israel_and_Cohen # [http://mp.cpgedupuydelome.fr/cours.php?id=3805&idPartie=34487] # - Gauss elimination is better.
# ======================================================================== # We are done print("You can run the file 'tests.py' to see examples of use of this module 'matrix.py'.") print("Testing every doctests in the module 'matrix'...") # Each function or method I wrote includes a doctest: import doctest doctest.testmod(verbose=True) doctest.testmod() print("\nMore details about doctest can be found on the Python documentation: \nhttps://docs.python.org/2/library/doctest.html") A = Matrix([[1, 2], [3, 4]])
# End of matrix.py |