Matrix Multiplication in pure Python

Question

I m trying to multiply two matrices together using pure Python  Input  X1 is a 3x3 and Xt is a 3x2    X1      1 0016  0 0  -16 0514            0 0  10000 0  -40000 0            -16 0514  -40000 0  160513 6437   Xt      1 0  1 0            0 0  0 25            0 0  0 0625     where Xt is the zip transpose of another matrix  Now here is the code   def matrixmult  A  B       C     0 for row in range len A    for col in range len B 0         for i in range len A            for j in range len B 0                 for k in range len B                    C i  j     A i  k  B k  j      return C   The error that python gives me is this       IndexError  list index out of range     Now I m not sure if Xt is recognised as an matrix and is still a list object  but technically this should work

User · Answer

def matrixmult  A  B       C     0 for row in range len A    for col in range len B 0         for i in range len A            for j in range len B 0                 for k in range len B                    C i  j     A i  k  B k  j      return C   at second line you should change  C     0 for row in range len B 0     for col in range len A

User · Answer

This is incorrect initialization  You interchanged row with col    C     0 for row in range len A    for col in range len B 0       Correct initialization would be  C     0 for col in range len B 0     for row in range len A      Also I would suggest using better naming conventions  Will help you a lot in debugging  For example   def matrixmult  A  B       rows A   len A      cols A   len A 0       rows B   len B      cols B   len B 0        if cols A    rows B        print  Cannot multiply the two matrices  Incorrect dimensions         return        Create the result matrix       Dimensions would be rows A x cols B     C     0 for row in range cols B   for col in range rows A       print C      for i in range rows A           for j in range cols B               for k in range cols A                   C i  j     A i  k    B k  j      return C   You can do a lot more  but you get the idea

User · Answer

All the below answers would return you the list Your need to convert it to matrix  def MATMUL X  Y       rows A   len X      cols A   len X 0       rows B   len Y      cols B   len Y 0        if cols A    rows B          print  Matrices are not compatible to Multiply  Check condition C1  R2          return        Create the result matrix       Dimensions would be rows A x cols B     C     0 for row in range cols B   for col in range rows A       print C      for i in range rows A           for j in range cols B               for k in range cols A                   C i  j     A i  k    B k  j       C   numpy matrix C  reshape len A  len B 0         return C

User · Answer

When I had to do some matrix arithmetic I defined a new class to help  Within such a class you can define magic methods like   add    or  in your use-case    matmul    allowing you to define x   a   b or a    b rather than matrixMult a b     matmul   was added in Python 3 5 per PEP 465   I have included some code which implements this below  I excluded the prohibitively long   init   method  which essentially creates a two-dimensional list self mat and a tuple self order according to what is passed to it   class Matrix      def   matmul   self  multiplier           if self order 1     multiplier order 0               raise ValueError  The multiplier was non-conformable under multiplication            return   sum a b for a b in zip srow mcol   for mcol in zip  multiplier mat   for srow in self mat       def   imatmul   self  multiplier           self mat   self   multiplier         return self mat      def   rmatmul   self  multiplicand           if multiplicand order 1     self order 0               raise ValueError  The multiplier was non-conformable under multiplication            return   sum a b for a b in zip mrow scol   for scol in zip  self mat   for mrow in multiplicand mat    Note      rmatmul   is used if b   a is called and b does not implement   matmul    e g  if I wanted to implement premultiplying by a 2D list    imatmul   is required for a    b to work correctly  If a matrix is non-conformable under multiplication it means that it cannot be multiplied  usually because it has more or less rows than there are columns in the multiplicand

User · Answer

Matrix Multiplication in pure python   def matmult m1 m2   r    m    for i in range len m1        for j in range len m2 0             sums 0         for k in range len m2                sums sums  m1 i  k  m2 k  j           r append sums      m append r      r    return m

User · Answer

One liner   def dot m1  m2       return            sum x   y for x  y in zip m1 r  m2 c   for m2 c in zip  m2   for m1 r in m1         Explanation   zip  m2  - gets a column from the second matrix  zip m1 r  m2 c  - creates tuple from m1 row and m2 column  sum      - sums multiplication row   col  Test   m1     1  2  3    4  5  6   m2     7  8    9  10    11  12   result   dot m1  m2  assert result      58  64    139  154

User · Answer

Here s some short and simple code for matrix vector routines in pure Python that I wrote many years ago      Basic Table  Matrix and Vector functions for Python 2 2    Author    Raymond Hettinger          Version    File MATFUNC PY  Ver 183  Date 12-Dec-2002 14 33 42   import operator  math  random NPRE  NPOST   0  0                      Disables pre and post condition checks  def iszero z    return abs z   lt   000001 def getreal z       try          return z real     except AttributeError          return z def getimag z       try          return z imag     except AttributeError          return 0 def getconj z       try          return z conjugate       except AttributeError          return z   separator           t     n     n---------- n     n            n     class Table list       dim   1     concat   list   add          A substitute for the overridden   add   method     def   getslice    self  i  j            return self   class    list   getslice   self i j        def   init    self  elems            list   init    self  elems           if len elems  and hasattr elems 0    dim    self dim   elems 0  dim   1     def   str    self            return separator self dim  join  map str  self        def map  self  op  rhs None               Apply a unary operator to every element in the matrix or a binary operator to corresponding         elements in two arrays   If the dimensions are different  broadcast the smaller dimension over         the larger  i e  match a scalar to every element in a vector or a vector to a matrix              if rhs is None                                                    Unary case             return self dim  1 and self   class    map op  self    or self   class     elem map op  for elem in self            elif not hasattr rhs  dim                                         List   Scalar op             return self   class     op e rhs  for e in self            elif self dim    rhs dim                                          Same level Vec   Vec or Matrix   Matrix             assert NPRE or len self     len rhs    Table operation requires len sizes to agree              return self   class    map op  self  rhs            elif self dim  lt  rhs dim                                           Vec   Matrix             return self   class     op self e  for e in rhs             return self   class     op e rhs  for e in self              Matrix   Vec     def   mul    self  rhs     return self map  operator mul  rhs       def   div    self  rhs     return self map  operator div  rhs       def   sub    self  rhs     return self map  operator sub  rhs       def   add    self  rhs     return self map  operator add  rhs       def   rmul    self  lhs     return self lhs     def   rdiv    self  lhs     return self  1 0 lhs      def   rsub    self  lhs     return - self-lhs      def   radd    self  lhs     return self lhs     def   abs    self    return self map  abs       def   neg    self    return self map  operator neg       def conjugate  self    return self map  getconj       def real  self    return self map  getreal        def imag  self    return self map  getimag       def flatten  self            if self dim    1  return self         return reduce  lambda cum  e  e flatten   concat cum   self           def prod  self     return reduce operator mul  self flatten    1 0      def sum  self     return reduce operator add  self flatten    0 0      def exists  self  predicate            for elem in self flatten                if predicate elem                   return 1         return 0     def forall  self  predicate            for elem in self flatten                if not predicate elem                   return 0         return 1     def   eq    self  rhs     return  self - rhs  forall  iszero    class Vec Table       def dot  self  otherVec     return reduce operator add  map operator mul  self  otherVec   0 0      def norm  self     return math sqrt abs  self dot self conjugate           def normalize  self     return self   self norm       def outer  self  otherVec     return Mat  otherVec x for x in self       def cross  self  otherVec             Compute a Vector or Cross Product with another vector          assert len self     len otherVec     3   Cross product only defined for 3-D vectors          u  v   self  otherVec         return Vec   u 1  v 2 -u 2  v 1   u 2  v 0 -u 0  v 2   u 0  v 1 -u 1  v 0         def house  self  index             Compute a Householder vector which zeroes all but the index element after a reflection          v   Vec  Table  0  index  concat self index      normalize           t   v index          sigma   1 0 - t  2         if sigma    0 0              t   v index    t lt  0 and t-1 0 or -sigma    t   1 0              v    t         return v  2 0   t  2    sigma   t  2      def polyval  self  x             Vec  6 3 4   polyval 5  evaluates to 6 x  2   3 x   4 at x 5          return reduce  lambda cum c  cum x c  self  0 0       def ratval  self  x             Vec  10 20 30 40 50   ratfit 5  evaluates to  10 x  2   20 x   30     40 x  2   50 x   1  at x 5           degree   len self    2         num  den   self  degree 1   self degree 1      1          return num polyval x    den polyval x   class Matrix Table         slots       size    rows    cols       def   init    self  elems             Form a matrix from a list of lists or a list of Vecs          Table   init    self  hasattr elems 0    dot   and elems or map Vec map tuple elems             self size   self rows  self cols   len elems   len elems 0       def tr  self             Tranpose elements so that Transposed i  j    Original j  i           return Mat zip  self       def star  self             Return the Hermetian adjoint so that Star i  j    Original j  i  conjugate            return self tr   conjugate       def diag  self             Return a vector composed of elements on the matrix diagonal          return Vec   self i  i  for i in range min self size          def trace  self    return self diag   sum       def mmul  self  other             Matrix multiply by another matrix or a column vector           if other dim  2  return Mat  map self mmul  other tr      tr           assert NPRE or self cols    len other          return Vec  map other dot  self        def augment  self  otherMat             Make a new matrix with the two original matrices laid side by side          assert self rows    otherMat rows   Size mismatch   s    s      self size    otherMat size           return Mat  map Table concat  self  otherMat        def qr  self  ROnly 0             QR decomposition using Householder reflections  Q R  self  Q tr   Q  I n   R upper triangular          R   self         m  n   R size         for i in range min m n                v  beta   R tr   i  house i              R -  v outer  R tr   mmul v  beta           for i in range 1 min n m    R i   i     0    i         R   Mat R  n           if ROnly  return R         Q   R tr   solve self tr    tr           Rt Qt   At    nn  nm    nm         self qr   lambda r 0  c  self   not r and c   self  and  Q R  or Matrix qr self r   Cache result         assert NPOST or m gt  n and Q size   m n  and isinstance R UpperTri  or m lt n and Q size   m m  and R size   m n          assert NPOST or Q mmul R   self and Q tr   mmul Q   eye min m n           return Q  R     def  solve  self  b               General matrices  incuding  are solved using the QR composition          For inconsistent cases  returns the least squares solution            Q  R   self qr           return R solve  Q tr   mmul b        def solve  self  b             Divide matrix into a column vector or matrix and iterate to improve the solution          if b dim  2  return Mat  map self solve  b tr      tr           assert NPRE or self rows    len b    Matrix row count  d must match vector length  d     self rows  len b           x   self  solve  b           diff   b - self mmul x          maxdiff   diff dot diff          for i in range 10               xnew   x   self  solve  diff               diffnew   b - self mmul xnew              maxdiffnew   diffnew dot diffnew              if maxdiffnew  gt   maxdiff   break             x  diff  maxdiff   xnew  diffnew  maxdiffnew              print  gt  gt  sys stderr  i 1  maxdiff         assert NPOST or self rows  self cols or self mmul x     b         return x     def rank  self     return Vec   not row forall iszero  for row in self qr ROnly 1     sum    class Square Matrix       def lu  self             Factor a square matrix into lower and upper triangular form such that L mmul U   A          n   self rows         L  U   eye n   Mat self             for i in range n               for j in range i 1 U rows                   assert U i  i     0 0   LU requires non-zero elements on the diagonal                  L j  i    m   1 0   U j  i    U i  i                  U j  -  U i    m         assert NPOST or isinstance L LowerTri  and isinstance U UpperTri  and L U  self         return L  U     def   pow    self  exp             Raise a square matrix to an integer power  i e  A  3 is the same as A mmul A mmul A            assert NPRE or exp  int exp  and exp gt 0   Matrix powers only defined for positive integers not  s    exp         if exp    1  return self         if exp amp 1  return self mmul self     exp-1           sqrme   self     exp 2          return sqrme mmul sqrme      def det  self     return self qr  ROnly 1   det       def inverse  self     return self solve  eye self rows        def hessenberg  self               Householder reduction to Hessenberg Form  zeroes below the diagonal          while keeping the same eigenvalues as self             for i in range self cols-2               v  beta   self tr   i  house i 1              self -  v outer  self tr   mmul v  beta               self -  self mmul v  outer v beta          return self     def eigs  self             Estimate principal eigenvalues using the QR with shifts method          origTrace  origDet   self trace    self det           self   self hessenberg           eigvals   Vec             for i in range self rows-1 0 -1               while not self i   i  forall iszero                   shift   eye i 1    self i  i                  q  r    self - shift  qr                   self   r mmul q    shift             eigvals append  self i  i                self   Mat   self r   i  for r in range i             eigvals append  self 0  0            assert NPOST or iszero   abs origDet  - abs eigvals prod       1000 0           assert NPOST or iszero  origTrace - eigvals sum             return Vec eigvals   class Triangular Square       def eigs  self     return self diag       def det  self     return self diag   prod    class UpperTri Triangular       def  solve  self  b             Solve an upper triangular matrix using backward substitution          x   Vec             for i in range self rows-1  -1  -1               assert NPRE or self i  i    Backsub requires non-zero elements on the diagonal              x insert 0   b i  - x dot self i  i 1       self i  i            return x  class LowerTri Triangular       def  solve  self  b             Solve a lower triangular matrix using forward substitution          x   Vec             for i in range self rows               assert NPRE or self i  i    Forward sub requires non-zero elements on the diagonal              x append   b i  - x dot self i   i      self i  i            return x  def Mat  elems         Factory function to create a new matrix       m  n   len elems   len elems 0       if m    n  return Matrix elems      if n  lt   1  return Square elems      for i in range 1  len elems            if not iszero  max map abs  elems i   i                   break     else  return UpperTri elems      for i in range 0  len elems -1           if not iszero  max map abs  elems i  i 1                    return Square elems      return LowerTri elems    def funToVec  tgtfun  low -1  high 1  steps 40  EqualSpacing 0           Compute x y points from evaluating a target function over an interval  low to high      at evenly spaces points or with Chebyshev abscissa spacing  default          if EqualSpacing          h    0 0 high-low  steps         xvec    low h 2 0 h i for i in range steps       else          scale  base    0 0 high-low  2 0   0 0 high low  2 0         xvec    base scale math cos   2 steps-1-2 i  math pi   2 steps   for i in range steps       yvec   map tgtfun  xvec      return Mat   xvec  yvec     def funfit   xvec  yvec   basisfuns         Solves design matrix for approximating to basis functions      return Mat   map form xvec  for form in basisfuns    tr   solve Vec yvec    def polyfit   xvec  yvec   degree 2         Solves Vandermonde design matrix for approximating polynomial coefficients      return Mat    x  n for n in range degree -1 -1   for x in xvec    solve Vec yvec    def ratfit   xvec  yvec   degree 2         Solves design matrix for approximating rational polynomial coefficients  a x  2   b x   c   d x  2   e x   1       return Mat   x  n for n in range degree -1 -1    -y x  n for n in range degree 0 -1   for x y in zip xvec yvec    solve Vec yvec    def genmat m  n  func       if not n  n m     return Mat    func i j  for i in range n   for j in range m      def zeroes m 1  n None        Zero matrix with side length m-by-m or m-by-n       return genmat m n  lambda i j  0   def eye m 1  n None        Identity matrix with side length m-by-m or m-by-n      return genmat m n  lambda i j  i  j   def hilb m 1  n None        Hilbert matrix with side length m-by-m or m-by-n   Elem i  j  1  i j 1       return genmat m n  lambda i j  1 0  i j 1 0    def rand m 1  n None        Random matrix with side length m-by-m or m-by-n      return genmat m n  lambda i j  random random     if   name         main         import cmath     a   Table  1 2j 2 3 4       b   Table  5 6 7 8       C   Table  a b       print  a b   a b     print  2 a   2 a     print  a 5 0   a 5 0     print  2 a 3 b   2 a 3 b     print  a C   a C     print  3 C   3 C     print  C b   C b     print  C sum     C sum       print  C map math cos    C map cmath cos      print  C conjugate     C conjugate       print  C real     C real        print zeroes 3      print eye 4      print hilb 3 5       C   Mat    1 2 3    4 5 1     7 8 9         print C mmul  C tr       n      print C    5    n      print C   C tr      n       A   C tr   augment  Mat   10 11 13    tr     tr       q  r   A qr       assert q mmul r     A     assert q tr   mmul q   eye 3      print  q  n   q    nr  n   r    nQ tr   amp Q  n   q tr   mmul q     nQ R n   q mmul r     n      b   Vec  50  100  220  321       x   A solve b      print  x      x     print  b      b     print  Ax     A mmul x       inv   C inverse       print   ninverse C  n   inv    nC   inv C   n   C mmul inv      assert C mmul inv     eye 3       points    xvec yvec    funToVec lambda x  math sin x  2 math cos  7 x  1   low 0  high 3  EqualSpacing 1      basis    lambda x  math sin x   lambda x  math exp x   lambda x  x  2      print  Func coeffs    funfit  points  basis       print  Poly coeffs    polyfit  points  degree 5       points    xvec yvec    funToVec lambda x  math sin x  2 math cos  7 x  1   low 0  high 3      print  Rational coeffs    ratfit  points        print polyfit   1 2 3 4    1 4 9 16    2       mtable   Vec  1 2 3   outer Vec  1 2        print mtable  mtable size      A   Mat    2 0 3    1 5 1    18 0 6         print  A       print A     print  eigs       print A eigs       print  Should be    Vec  11 6158  5 0000  -3 6158       print  det A       print A det        c   Mat    1 2 30   4 5 10   10 80 9           Failed example from Konrad Hinsen     print  C  n   c     print c eigs       print  Should be    Vec  -8 9554  43 2497  -19 2943        A   Mat    1 2 3 4    4 5 6 7    2 1 5 0    4 2 1 0           Kincaid and Cheney p 326     print  A  n   A     print A eigs       print  Should be    Vec  3 5736  0 1765  11 1055  -3 8556        A   rand 3      q r   A qr       s t   A qr       print q is s                  Test caching     print r is t     A 1  1    1 1                 Invalidate the cache     u v   A qr       print q is u                  Verify old result not used     print r is v     print u mmul v     A          Verify new result      print  Test qr on 3x5 matrix      a   rand 3 5      q r   a qr       print q mmul r     a     print q tr   mmul q     eye 3

User · Answer

If you really don t want to use numpy you can do something like this   def matmult a b       zip b   zip  b        uncomment next line if python 3          zip b   list zip b      return   sum ele a ele b for ele a  ele b in zip row a  col b                 for col b in zip b  for row a in a   x     1 2 3   4 5 6   7 8 9   10 11 12   y     1 2   1 2   3 4    import numpy as np   I want to check my solution with numpy  mx   np matrix x  my   np matrix y           Result    gt  gt  gt  matmult x y    12  18    27  42    42  66    57  90    gt  gt  gt  mx   my matrix   12  18            27  42            42  66            57  90

User · Answer

The fault occurs here   C i  j   A i  k  B k  j    It crashes when k 2  This is because the tuple A i  has only 2 values  and therefore you can only call it up to A i  1  before it errors   EDIT  Listen to Gerard s answer too  your C is wrong  It should be C   0 for row in range len A    for col in range len A 0       Just a tip  you could replace the first loop with a multiplication  so it would be C   0  len A  for col in range len A 0

User · Answer

m input  row   n input  col   X    for i in range  m       m1        for j in range  n           m1 append input  num        X append m1  Y    for i in range  m       n1        for j in range  n           n1 append input  num        Y append n1      result is 3x3 result     0 0 0             0 0 0             0 0 0      iterate through rows of X for i in range len X         iterate through columns of Y    for j in range len Y 0              iterate through rows of Y        for k in range len Y               result i  j     X i  k    Y k  j   for r in result     print r

User · Answer

The shape of your matrix C is wrong  it s the transpose of what you actually want it to be   But I agree with ulmangt  the Right Thing is almost certainly to use numpy  really

[python] Matrix Multiplication in pure Python?

Examples related to python

Examples related to matrix-multiplication