#! /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