Source code for matrix

#! /usr/bin/env python
# -*- coding: utf-8 -*-
""" Complete solution for the CS101 Programming Project about matrices.

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>`_.
"""

from __future__ import division, print_function, absolute_import  # Python 2/3 compatibility

from math import factorial as _factorial  # For exp(A)


# Experimental: computation with decimal numbers to improve decimal precision
from decimal import Decimal as _Decimal


[docs]class Decimal(_Decimal): """ 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. """ def __str__(self, *args, **kwargs): if int(self) == self: return "{}".format(int(self)) else: return _Decimal.__str__(self) __repr__ = __str__
# Experimental: computation with fraction to be exact and not numerically approximative ! from fractions import Fraction as _Fraction
[docs]class Fraction(_Fraction): """ Extended :class:`fractions.Fraction` class to improve the ``str`` and ``repr`` methods. If the denominator is 1, print it as an integer. """ def __str__(self, *args, **kwargs): if self.denominator == 1: return "{}".format(self.numerator) else: return "{}/{}".format(self.numerator, self.denominator) __repr__ = __str__
# ========================================================================
[docs]class Matrix(object): """ 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. """
[docs] def __init__(self, listrows): """ 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]] """ try: self._n = len(listrows) if self._n == 0: self._m = 0 else: self._m = len(listrows[0]) assert all(self.m == len(listrows[i]) for i in range(self._n)) # Now we do get a fresh copy of that list. #: self.listrows is the list of rows for self self.listrows = [[listrows[i][j] for j in range(self._m)] for i in range(self._n)] # FIXME We should forbid modifying these attributes from outside the class self._i0, self._j0 = 0, 0 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 @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 """ return self._n rows = n # This decorator @property makes this method an attributes # cf. https://docs.python.org/2/library/functions.html#property @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 """ return self._m cols = m # ==================================================================== # Methods for reading and accessing
[docs] def __getitem__(self, ij): """ ``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, j = ij if isinstance(i, int): if isinstance(j, int): return self.listrows[i][j] elif isinstance(j, slice): # i is an integer, j is a slice object return Matrix([self.listrows[i][j]]) elif isinstance(i, slice): if isinstance(j, int): # i is a slice object, j is an integer return Matrix([[x[j]] for x in self.listrows[i]]) elif isinstance(j, slice): # i and j are a slice objects return Matrix([x[j] for x in self.listrows[i]]) # 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)))
[docs] def __setitem__(self, ij, value): """ ``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]] """ i, j = ij fail = False if isinstance(value, Matrix): value = value.listrows # just the list of rows if isinstance(i, int): if isinstance(j, int): self.listrows[i][j] = value # This is the simple case, the one we use the most elif isinstance(j, slice): # i is an integer, j is a slice object j_value = 0 # for j0 in _slice_to_range(j): for j0 in range(*j.indices(self.cols)): try: self.listrows[i][j0] = value[0][j_value] # sub-column except Exception: try: self.listrows[i][j0] = value[j_value] # list except Exception: self.listrows[i][j0] = value j_value += 1 # End for loop j0 else: fail = True elif isinstance(i, slice): if isinstance(j, int): # i is a slice object, j is an integer i_value = 0 # for i0 in _slice_to_range(i): for i0 in range(*i.indices(self.rows)): try: self.listrows[i0][j] = value[i_value][0] # sub-row except Exception: try: self.listrows[i0][j] = value[i_value] # list except Exception: self.listrows[i0][j] = value i_value += 1 # 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 if fail: # 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
[docs] def row(self, i): """ ``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]] """ return Matrix([[self[i, j] for j in range(self.m)]])
[docs] def col(self, j): """ ``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]] """ return Matrix([[self[i, j]] for i in range(self.n)])
# ==================================================================== # 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)
[docs] def copy(self): """ ``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]] """ return Matrix(self.listrows)
# ==================================================================== # Length and shape
[docs] def __len__(self): """ ``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 """ return self.n * self.m
# This decorator @property makes this method an attributes # cf. https://docs.python.org/2/library/functions.html#property @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) """ return (self.n, self.m) # ==================================================================== # Transposition
[docs] def transpose(self): """ ``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 """ return Matrix([[self[j, i] for j in range(self.n)] for i in range(self.m)])
@property 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 """ return self.transpose() # ==================================================================== # Methods for pretty-printing
[docs] def __str__(self): """ ``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]]' """ try: return str(self.map(str).listrows).replace("'", "") except Exception: str(self.listrows)
[docs] def __repr__(self): """ ``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 ==
[docs] def __eq__(self, B): 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 """ 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
[docs] def almosteq(self, B, epsilon=1e-10): 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 """ try: if self.n == B.n and self.m == B.m: return all(abs(self[i, j] - B[i, j]) < epsilon for j in range(self.m) for i in range(self.n)) else: return False except Exception: return False
# Comparing <
[docs] def __lt__(self, B): 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)
[docs] def __add__(self, B): 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]] """ if isinstance(B, Matrix): # Sum of two matrices assert self.n == B.n and self.m == B.m return Matrix([[self[i, j] + B[i, j] for j in range(self.m)] for i in range(self.n)]) else: # Sum of matrix A and a number B return Matrix([[self[i, j] + B for j in range(self.m)] for i in range(self.n)])
[docs] def __radd__(self, 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]] """ if isinstance(B, Matrix): # 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) return Matrix([[B + self[i, j] for j in range(self.m)] for i in range(self.n)])
# ==================================================================== # Substraction (left and right)
[docs] def __sub__(self, B): 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]] """ if isinstance(B, Matrix): # Sum of two matrices assert self.n == B.n and self.m == B.m return Matrix([[self[i, j] - B[i, j] for j in range(self.m)] for i in range(self.n)]) else: # Sum of matrix A and a number B return Matrix([[self[i, j] - B for j in range(self.m)] for i in range(self.n)])
[docs] def __neg__(self): 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 """ return Matrix([[-self[i, j] for j in range(self.m)] for i in range(self.n)])
[docs] def __pos__(self): 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 """ return Matrix([[+self[i, j] for j in range(self.m)] for i in range(self.n)])
[docs] def __rsub__(self, B): 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 """ if isinstance(B, Matrix): # 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) return Matrix([[B - self[i, j] for j in range(self.m)] for i in range(self.n)])
# ==================================================================== # Product (left and right)
[docs] def __mul__(self, B): 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]] """ if isinstance(B, Matrix): # Product of two matrices assert self.m == B.n return Matrix([[sum(self[i, k] * B[k, j] for k in range(self.m)) for j in range(B.m)] for i in range(self.n)]) else: # Product 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)])
[docs] def __rmul__(self, B): 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]] """ if isinstance(B, Matrix): # 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) return Matrix([[B * self[i, j] for j in range(self.m)] for i in range(self.n)])
[docs] def multiply_elementwise(self, B): 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]] """ if not isinstance(B, Matrix): raise ValueError("A.multiply_elementwise(B): B has to be a Matrix object.") else: assert self.shape == B.shape return Matrix([[self[i, j] * B[i, j] for j in range(self.m)] for i in range(self.n)])
# ==================================================================== # Division (left and right)
[docs] def __div__(self, B): 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. 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] / float(B) for j in range(self.m)] for i in range(self.n)]) return Matrix([[self[i, j] / B for j in range(self.m)] for i in range(self.n)])
__truediv__ = __div__
[docs] def __floordiv__(self, B): 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)])
[docs] def __mod__(self, b): 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]] """ if isinstance(b, Matrix): # 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: return Matrix([[self[i, j] % b for j in range(self.m)] for i in range(self.n)])
[docs] def __rdiv__(self, B): 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)])
__rtruediv__ = __rdiv__
[docs] def __rfloordiv__(self, B): 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
[docs] def __pow__(self, k): 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 """ if not isinstance(k, int): raise ValueError("A ** k: k should be an integer (here k = {}).".format(k)) elif k < 0: # A ^ k = (A ^ (-1)) ^ (-k) return (self.inv()) ** (abs(k)) elif k == 0: return eye(self.n) # This is a convention : a ** 0 = 1_n,n elif k == 1: return self.copy() # It is really import to get a copy ! elif k == 2: # Useless in fact return self * self elif k % 2 == 0: P = self * self return P ** int(k // 2) # A^(2k) = ((A**2) ^ k) elif k % 2 == 1: P = self * self return self * (P ** int((k - 1) // 2)) # A^(2k+1) = A * ((A**2) ^ k) # 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))
[docs] def exp(self, limit=30): 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 """ assert limit > 0 if not self.is_square: raise ValueError("A.exp() is only possible if A is a square matrix.") e = eye(self.n) for k in range(1, limit): e += (self ** k) * (1.0 / _factorial(k)) return e
[docs] def inv(self): """ ``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 """ if not self.is_square: raise ValueError("A.inv() is only possible if A is a square matrix.") else: try: _, inverse = self.gauss_jordan(inv=True) return inverse except Exception: raise ValueError("A.inv() on a singular matrix (ie. non inversible).")
# ==================================================================== # Gauss elimination process (to get a row echelon form)
[docs] def gauss(self, det=False, verb=False, mode=None, maxpivot=False): 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. c = self.copy() n, m = c.n, c.m currentdet = 1 if mode is None: # Trying to compute the mode ourself if all(isinstance(x, Fraction) for x in self): mode = 'f' elif all(isinstance(x, Decimal) for x in self): 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() for k in range(min(n, m)): if verb: print("\nTrying to find the {}-th pivot:".format(k)) print("With these indeces:", list(range(k, m))) print("And that array:", [abs(c[k, j]) for j in range(m)]) if maxpivot: i_max = _argmax(list(range(k, m)), [abs(c[k, j]) for j in range(m)]) else: i_max = k if c[k, i_max] == 0: for possible_i_max in range(k, m): if c[k, possible_i_max] != 0: i_max = possible_i_max break # We found the first i_max such that c[k, i_max] is not zero, or k if none are good if verb: # assert c[k, i_max] == max(abs(c[k, j]) for j in range(m)) print("_argmax has given i_max = {}, and c[k, i_max] = {} (with k = {}).".format(i_max, c[k, i_max], k)) if c[k, i_max] == 0: currentdet = 0 # XXX Do we really have a singular matrix already ? # raise ValueError("A.gauss_elimination() called on a singular matrix.") if (not det) and verb: print("WARNING: A.gauss() might have been called on a singular matrix. FIXME remove these warnings") # return (c, 0) if det else c # XXX not yet ! # determinant is 0, that is sure at least else: if verb: print("c.col(i_max) is:", c.col(i_max)) print("c.col(k) is:", c.col(k)) print("i_max = {}, k = {}.".format(i_max, k)) # c.swap_rows(i_max, k) c.swap_cols(i_max, k) # XXX Shouldn't we swap rows instead? I think not # Swapping two rows multiplies the determinant by -1 if i_max != k: if verb: print("We swapped two different lines ({} and {}), the determinant will be multiplied by -1.".format(i_max, k)) currentdet *= -1 if k >= (min(n, m) - 1): if verb: 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)) break # break the for loop RIGHT NOW if verb: print("Gauss Elimination: using the {k}th line (L_{k} = {l}).\n We use {pivot} as a pivot.".format(k=k, l=c.row(k), pivot=c[k, k])) # Do for all lines below pivot: for i in range(k + 1, n): if mode == 'f': gamma = Fraction(c[i, k], c[k, k]) elif mode == 'd': gamma = Decimal(c[i, k]) / Decimal(c[k, k]) else: # gamma = float(c[i, k]) / float(c[k, k]) gamma = c[i, k] / c[k, k] if verb: print(" Operation L_{i}' <-- L_{i} - gamma * L_{k}".format(i=i, k=k)) print(" with gamma =", gamma) print(" with old L_{i} = {l}".format(i=i, l=c.row(i))) # Do for all remaining elements in current line: for j in range(k + 1, m): c[i, j] -= gamma * c[k, j] # We convert to integer if possible, it is prettier :) # FIXME isn't it a cause of rounding mistake? if int(c[i, j]) == c[i, j]: c[i, j] = int(c[i, j]) # Fill lower triangular matrix with zeros (because gamma is chosen like that): c[i, k] = 0 if verb: print(" with new L_{i}' = {l}".format(i=i, l=c.row(i))) if det: # Product of the (-1)**(nb of swaps) * diagonal elements currentdet *= _prod(c[i, i] for i in range(min(n, m))) return c, currentdet else: return c
# End of gauss() # ==================================================================== # Gauss-Jordan elimination process (to get a reduced row echelon form)
[docs] def gauss_jordan(self, inv=False, verb=False, mode=None, maxpivot=False): """ ``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. c = self.copy() if mode is None: # Trying to compute the mode ourself if all(isinstance(x, Fraction) for x in self): mode = 'f' elif all(isinstance(x, Decimal) for x in self): mode = 'd' elif mode == 'f': try: c = self.map(Fraction) 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() if inv: if c.n != c.m: raise ValueError("A.gauss_jordan(inv=True) is only possible if A is a square matrix.") cinv = eye(c.n) # 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) n, m = c.n, c.m r = -1 for j in range(m): if verb: 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)]) if maxpivot: k = _argmax(list(range(r + 1, n)), [abs(c[i, j]) for i in range(n)]) else: k = r + 1 if c[k, j] == 0: 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 if verb: print("For r = {}, c[k, j] is the pivot (k = {}, j = {}), equals to {}.".format(r + 1, k, j, c[k, j])) if inv and c[k, j] == 0: raise ValueError("A.gauss_jordan() called on a singular matrix.") if c[k, j] != 0: pivot = c[k, j] r += 1 if verb: 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 for jjjj in range(m): # c[k, jjjj] /= pivot if mode == 'f': c[k, jjjj] = Fraction(c[k, jjjj], pivot) elif mode == 'd': c[k, jjjj] = Decimal(c[k, jjjj] / pivot) else: c[k, jjjj] /= pivot if inv: if verb: print("INFO: Same linear operation is applied to cinv: cinv[k, :] /= pivot") # cinv[k, :] /= pivot for jjjj in range(m): # cinv[k, jjjj] /= pivot if mode == 'f': cinv[k, jjjj] = Fraction(cinv[k, jjjj], pivot) elif mode == 'd': cinv[k, jjjj] = Decimal(cinv[k, jjjj] / pivot) else: cinv[k, jjjj] /= pivot if k != r: 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, :]) for i in range(n): if verb: print("For i = {}.".format(i)) if i != r and c[i, j] != 0: cij = c[i, j] if verb: 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)) c[i, :] -= c[r, :] * cij if inv: if verb: print("INFO: Same linear operation is applied to cinv: cinv[i, :] -= cinv[r, :] * c[i, j]") cinv[i, :] -= cinv[r, :] * cij # R_i <-- R_i - c[i,j] * R_r if verb: print("After: R_{i} =", c[i, :]) else: if verb: print("No modification here (i = r or c[i, j] = 0).") # Done if verb: print("Done ! c =", c) # c = c.map(lambda x: int(x) if int(x) == x else x) # Pretty priting only if inv: if verb: print("Done ! cinv =", cinv) return c, cinv else: return c
# ==================================================================== # Applications of the Gauss elimination process @property def rank(self, verb=False): """ ``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 """ c = self.gauss(verb=verb, det=False) return sum(c[i, i] != 0 for i in range(min(self.n, self.m))) @property def det(self, verb=False): """ ``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 """ _, d = self.gauss(det=True, verb=verb) if int(d) == d: d = int(d) return d # ==================================================================== # Extra methods (not required in the project)
[docs] def count(self, value): """ ``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 """ return sum(self[i, j] == value for j in range(self.m) for i in range(self.n))
[docs] def __contains__(self, value): """ ``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 """ return any(self[i, j] == value for j in range(self.m) for i in range(self.n))
[docs] def map(self, f, *args, **kwargs): """ 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]] """ return Matrix([[f(self[i, j], *args, **kwargs) for j in range(self.m)] for i in range(self.n)])
[docs] def round(self, ndigits=8): """ ``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
[docs] def __iter__(self): """ ``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] """ for i in range(self.n): for j in range(self.m): yield self[i, j]
# next and __next__ for iterating over the values of our matrix
[docs] def __next__(self): """ For Python 3 compatibility.""" return self.next()
[docs] def next(self): """ 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 @property 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]] """ return self.map(lambda z: z.real if isinstance(z, complex) else z) @property 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]] """ return self.map(lambda z: z.imag if isinstance(z, complex) else z)
[docs] def conjugate(self): """ Conjugate part of the matrix ``A``, coefficient wise. >>> A = Matrix([[-1j, -2j], [-2j, -1j]]) >>> A.conjugate() [[1j, 2j], [2j, 1j]] """ return self.map(lambda z: z.conjugate() if isinstance(z, complex) else z)
# ==================================================================== # Dot product and norm
[docs] def dot(self, v): 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. """ if isinstance(v, Matrix): if self.m == v.n and v.m == 1: return self * v 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 try: vector_v = Matrix([[x] for x in list(v)]) except Exception: raise ValueError(("A.dot(v): impossible to convert v = {} to a column vector.".format(v))) return self.dot(vector_v)
[docs] def norm(self, p=2): 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 """ if p == 'inf': return max(abs(x) for x in self) else: result = sum(abs(x)**p for x in self) ** (1.0 / p) return int(result) if int(result) == result else result
[docs] def normalized(self, fnorm=None, *args, **kwargs): """ ``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
[docs] def __abs__(self): """ ``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]] """ return self.map(abs)
abs = __abs__ # ==================================================================== # Trace and other special values
[docs] def trace(self): 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 """ return sum(self[i, i] for i in range(min(self.n, self.m)))
# 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 @property 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 """ return self.n == self.m @property 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 """ return self.n == self.m and all(self[i, j] == self[j, i] for i in range(self.n) for j in range(self.m)) @property 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 """ return self.n == self.m and all(self[i, j] == -self[j, i] for i in range(self.n) for j in range(self.m)) @property 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 """ return all(self[i, j] == 0 for i in range(self.n) for j in range(self.m) if i != j) @property 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 f(z): return z.conjugate() if isinstance(z, complex) else z return self.n == self.m and all(self[i, j] == f(self[j, i]) for i in range(self.n) for j in range(self.m)) @property 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 """ return all(self[i, j] == 0 for i in range(self.n) for j in range(i + 1, self.m)) @property 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 """ return all(self[i, j] == 0 for i in range(1, self.n) for j in range(i)) @property 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 """ return all(self[i, j] == 0 for j in range(self.m) for i in range(self.n)) @property 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*
[docs] def swap_cols(self, j1, j2): """ ``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]] """ for i in range(self.rows): self[i, j1], self[i, j2] = self[i, j2], self[i, j1]
[docs] def swap_rows(self, i1, i2): """ ``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]] """ for j in range(self.cols): self[i1, j], self[i2, j] = self[i2, j], self[i1, j]
# ======================================================================== # Adjugate matrix (https://en.wikipedia.org/wiki/Adjugate_matrix)
[docs] def minor(self, i, j): 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 """ return det(Matrix([[self[i0, j0] for j0 in range(self.m) if j0 != j] for i0 in range(self.n) if i0 != i]))
[docs] def cofactor(self, i, j): 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 """ return (-1)**(i + j) * self.minor(i, j)
[docs] def adjugate(self): 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 """ return mat_from_f(self.cofactor, self.n)
[docs] def type(self): """ ``A.type()`` returns the matrix of types of coefficients of ``A``. """ return self.map(type)
# End of that class Matrix # ======================================================================== # ======================================================================== # Utility functions
[docs]def ones(n, m=None): """ ``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]] """ if not m: m = n if isinstance(n, tuple): n, m = n return Matrix([[1 for _ in range(m)] for _ in range(n)])
[docs]def zeros(n, m=None): """ ``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]] """ if not m: m = n if isinstance(n, tuple): n, m = n return Matrix([[0 for _ in range(m)] for _ in range(n)])
[docs]def eye(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 """ return Matrix([[1 if i == j else 0 for j in range(n)] for i in range(n)])
[docs]def diag(d, n=None): """ ``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]] """ if n: return diag([d] * n) else: n = len(d) return Matrix([[d[i] if i == j else 0 for j in range(n)] for i in range(n)])
[docs]def mat_from_f(f, n, m=None, *args, **kwargs): """ ``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``). """ if not m: m = n return Matrix([[f(i, j, *args, **kwargs) for j in range(m)] for i in range(n)])
# ======================================================================== # Functions that are just wrappers around methods
[docs]def det(A): 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 """ return A.det
[docs]def rank(A): r""" ``rank(A)`` <-> ``A.rank`` computes the rank of ``A`` (in :math:`\mathcal{O}(n^3)`). >>> rank(eye(2)) 2 """ return A.rank
[docs]def gauss(A, *args, **kwargs): r""" ``gauss(A)`` <-> ``A.gauss()`` applies the Gauss elimination process on ``A`` (in :math:`\mathcal{O}(n^3)`). """ return A.gauss_elimination(*args, **kwargs)
[docs]def gauss_jordan(A, *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)
[docs]def inv(A): r""" ``inv(A)`` <-> ``A.inv()`` **tries** to compute the inverse of ``A`` (in :math:`\mathcal{O}(n^3)`). >>> inv(eye(2)) == eye(2) True """ return A.inv()
[docs]def exp(A, *args, **kwargs): 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 """ return A.exp(*args, **kwargs)
# ======================================================================== # The LU decomposition
[docs]def PLUdecomposition(A, mode=None): 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). """ assert A.is_square U = A.copy() n, m = U.n, U.m P, L = eye(n), eye(n) if mode is None: # Trying to compute the mode ourself if all(isinstance(x, Fraction) for x in U): mode = 'f' 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 for k in range(n): # Find the pivot on the k-th row i_max = k # Is the pivot U[k, k] is zero, we find a possible better pivot if U[k, i_max] == 0: i_max = _argmax(list(range(k, m)), [abs(U[k, j]) for j in range(m)]) # Is the pivot U[k, i_max] is still non zero, we cannot do anything, because the matrix is singular. if U[k, i_max] == 0: raise ValueError("PLUdecomposition() has been called on a singular matrix.") else: # U.swap_cols(i_max, k) U.swap_rows(i_max, k) # The matrix P will keep track of the permutations performed during the Gaussian Elimination process: # P.swap_cols(i_max, k) P.swap_rows(i_max, k) if k >= (min(n, m) - 1): break # break the for loop RIGHT NOW # Do for all lines/rows below pivot: for i in range(k + 1, n): if mode == 'f': gamma = Fraction(U[i, k], U[k, k]) 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: for j in range(k + 1, m): # Add - gamma times row k to row i of U U[i, j] -= gamma * U[k, j] # 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): if mode == 'f': U[i, k] = Fraction(0) 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. L[i, k] = gamma # 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. for j in range(min(i_max, k) + 1, m): L[i_max, j], L[k, j] = L[k, j], L[i_max, j] # L[j, i_max], L[j, k] = L[j, k], L[j, i_max] return P, L, U
# End of PLUdecomposition() # ======================================================================== # Other functions
[docs]def norm(A, p=2, *args, **kwargs): """ ``norm(A, p)`` <-> ``A.norm(p)`` computes the p-norm of ``A`` (default is ``p = 2``).""" return A.norm(p, *args, **kwargs)
[docs]def trace(A, *args, **kwargs): """ ``trace(A)`` <-> ``A.trace()`` computes the trace of ``A``.""" return A.trace(*args, **kwargs)
# ======================================================================== # Random matrix from random import randint as _randint, uniform as _uniform
[docs]def rand_matrix(n=1, m=1, k=10): """ ``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]] """ assert k > 0 return Matrix([[_randint(-k, k) for _ in range(m)] for _ in range(n)])
[docs]def rand_matrix_float(n=1, m=1, k=10): """ ``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]] """ assert k > 0 return Matrix([[_uniform(-k, k) for _ in range(m)] for _ in range(n)])
# ======================================================================== # 2 auxiliary functions used by the Gauss elimination process
[docs]def _argmax(indexes, array): """ Compute the index ``i`` in ``indexes`` such that the ``array[i]`` is the bigger.""" besti = None if not indexes: raise ValueError("argmax() arg indexes is a non-valid sequence.") # bestvalue = array[indexes[0]] bestvalue = float('-inf') # Comparison with None fails in Python 3 for i in indexes: if array[i] >= bestvalue: # WARNING Accessing array[i] does not check if i is a good index or not besti = i bestvalue = array[i] if besti is None: raise ValueError("argmax() arg is a non-valid sequence.") return besti
[docs]def _prod(iterator): """ Compute the product of the values in the iterator ``iterator``. Empty product is 1.""" p = 1 for x in iterator: p *= x return p
# 2 auxiliary functions for implementing the generalized __setitem__ method
[docs]def _ifnone(a, b): """ ``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
[docs]def _slice_to_range(sliceobject): """ 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
[docs]def innerproduct(vx, vy): 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))
[docs]def norm_square(u): 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
[docs]def norm2(u): 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
[docs]def vect_const_multi(vx, c): """ 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]
[docs]def proj(u, v): 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)
[docs]def gram_schmidt(V, normalized=False): 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)
[docs]def minor(A, i, j): 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``). """ return A.minor(i, j)
[docs]def cofactor(A, i, j): 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``). """ return A.cofactor(i, j)
[docs]def adjugate(A): 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)`. """ return A.adjugate()
# ======================================================================== # 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 if __name__ == '__main__': 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