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.

Realgar

@ulmangt: "using pure python". He/she wants to do it without downloadable modules, probably for the challenge. – Candiecandied

@ulmangt, not all implementations of Python can use numpy/scipy – Ladida

Yeah...the challenge...Thanks @beary605. – Realgar

By the way, from Python 3.5 a special operator '@' can be used for matrix multiplication (such as X @ W + b). legacy.python.org/dev/peps/pep-0465 – Tramontane

@agf The homework tag is deprecated, see meta.stackexchange.com/questions/147100/… – Angling

@Angling But it wasn't when I posted that. – Concinnous

@Concinnous ... True ... – Angling

@ulmangt, not all installations of Python are on machines where extensions like numpy that are written C can be built. – Thereupon

Answer 1

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:

>>> matmult(x,y)
[[12, 18], [27, 42], [42, 66], [57, 90]]
>>> mx * my
matrix([[12, 18],
        [27, 42],
        [42, 66],
        [57, 90]])

Tryparsamide

You can easily improve this by only computing zip(*b) once – Ladida

@gnibbler, that's a good point, thank you. I edited the code to reflect your suggestion. – Tryparsamide

I think you can replace the whole summation with sum(map(mul, row_a, col_b)) (after doing from operator import mul, if you don't want to import, sum(map(lambda x, y: x * y, row_a, col_b)). – Plutonian

Answer 2

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...

Feinleib

Answer 3

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) < .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 < 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<=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>=n and Q.size==(m,n) and isinstance(R,UpperTri) or m<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 >= maxdiff:  break
            x, diff, maxdiff = xnew, diffnew, maxdiffnew
            #print >> 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>0, 'Matrix powers only defined for positive integers not %s' % exp
        if exp == 1: return self
        if exp&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 <= 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()&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)

Wanting

Answer 4

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]]

Emergency

Answer 5

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

Doughnut

For those not always working with the latest versions of Python: this matrix multiplication operator was added in Python 3.5. – Deeplaid

For details on the Python 3.5 addition, see "PEP 465 -- A dedicated infix operator for matrix multiplication" python.org/dev/peps/pep-0465 – Fredkin

Answer 6

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

Upheld

Answer 7

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]))]

Candiecandied

True if matrixMult(Xt,X1) is evaluated – Fatma

Answer 8

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.)

Gaidano

Answer 9

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

Oxalate

Answer 10

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))]

Bedew

Answer 11

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)

Counterblast

Answer 12

Had to delete the first post because I noticed an error but this should work fine for multiplying two matrixes.

from numpy import*


A = matrix([[1.0016, 0.0, -16.0514],[0.0, 10000.0, -40000.0],[-16.0514, -40000.0, 160513.6437]])
B = matrix([(1.0, 1.0),(0.0, 0.25),(0.0, 0.0625)])

n,m = shape(B)
C = ones((n,m*n))
index=0
matrixA=ones(shape(A))

for i in range(n):
    if (i==1):
        for y in range(n):
            for z in range(n):
                matrixA[y,z]=matrixA[y,z]*A[y,z]
    elif(i>1):
        matrixA=matrixA*A
    if (i==0):
        result=B
    else:
        result = matrixA * B
    for x in range (n):
        for y in range(m):
            C[x,index] = result[x,y]
            index+=1
        index=index-m
    index=index+m


print(C)

Benoni