transpuesta - ¿Multiplicación de matrices en python?
suma de matrices en python (10)
Aquí hay algunos códigos cortos y simples para rutinas de matrices / vectores en Python puro que escribí hace muchos años:
''''''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)
Estoy tratando de multiplicar dos matrices usando python puro. Entrada (X1 es un 3x3 y Xt es un 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)]
donde Xt es la transposición zip de otra matriz. Ahora aquí está el código:
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
El error que Python me da es este: IndexError: índice de lista fuera de rango. Ahora no estoy seguro de si Xt se reconoce como una matriz y sigue siendo un objeto de lista, pero técnicamente esto debería funcionar.
Esta es una inicialización incorrecta. ¡Intercambiaste fila con col!
C = [[0 for row in range(len(A))] for col in range(len(B[0]))]
La inicialización correcta sería
C = [[0 for col in range(len(B[0]))] for row in range(len(A))]
También sugeriría usar mejores convenciones de nomenclatura. Te ayudará mucho en la depuración. Por ejemplo:
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
Puedes hacer mucho más, pero tienes la idea ...
La culpa se produce aquí:
C[i][j]+=A[i][k]*B[k][j]
Se bloquea cuando k = 2. Esto se debe a que la tupla A[i] tiene solo 2 valores y, por lo tanto, solo puede llamarla a A [i] [1] antes de que se produzca un error.
EDITAR: escuchar la respuesta de Gerard también, su C está mal. Debe ser C=[[0 for row in range(len(A))] for col in range(len(A[0]))] .
Solo una sugerencia: podría reemplazar el primer bucle con una multiplicación, por lo que sería C=[[0]*len(A) for col in range(len(A[0]))]
La forma de tu matriz C es incorrecta; Es la transposición de lo que realmente quieres que sea. (Pero estoy de acuerdo con ulmangt: lo correcto es, en verdad, usar NUMPY, en realidad).
Multiplicación de matrices en pitón puro.
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
Si realmente no quieres usar numpy puedes hacer algo como esto:
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)
Resultado:
>>> matmult(x,y)
[[12, 18], [27, 42], [42, 66], [57, 90]]
>>> mx * my
matrix([[12, 18],
[27, 42],
[42, 66],
[57, 90]])
Tarde a la fiesta, pero cuando tuve que hacer una aritmética matricial definí una nueva clase para ayudar. Dentro de dicha clase, puede definir métodos como __add__ , o, en su caso de uso, __matmul__ , lo que le permite hacer a * b o a *= b lugar de matrixMult(a,b) .
He incluido algún código que implementa esto a continuación ( __init__ método prohibitivamente largo __init__ , que esencialmente crea una lista bidimensional de self.mat y una tupla self.order acuerdo con lo que se le pasa)
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]
Nota:
- Se requiere
__rmatmul__para quea @ byb @ afuncionen correctamente; - Se requiere
__imatmul__para quea @= bfuncione correctamente; - Si una matriz no es conforme con la multiplicación, significa que no se puede multiplicar, generalmente porque tiene más o menos filas que columnas en el multiplicando
EDITAR: Me acabo de dar cuenta de que __rmatmul__ no es obligatorio ya que solo se usa para evaluar a @ b si a no tiene el método __matmul__ . Como solo estoy multiplicando las matrices por otras instancias de matrix , tendré un método __matmul__ , pero si alguna vez lo edito para que las operaciones funcionen, digamos, con matrices 2D, tendré que volver a agregar __rmatmul__ para que __rmatmul__ la list @ matrix así como la matrix @ list .
Todas las respuestas a continuación le devolverían la lista. Debe convertirla a matriz.
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
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
en segunda linea deberias cambiar
C = [[0 for row in range(len(B[0]))] for col in range(len(A))]
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)