blocks|key|2885013|text|你的矩阵C的形状是错误的；它是你真正想要的转置。(但我同意ulmangt:几乎可以肯定的是，使用numpy是正确的，真的。)|type|unstyled|depth|inlineStyleRanges|offset|length|style|CODE|entityRanges|data|2885014|entityMap^0|4|1|0^^$0|@$1|2|3|4|5|6|7|H|8|@$9|I|A|J|B|C]]|D|@]|E|$]]|$1|F|3|-4|5|6|7|K|8|@]|D|@]|E|$]]]|G|$]]

The shape of your matrix <code>C</code> 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.)

blocks|key|2885023|text|故障发生在以下位置：|type|unstyled|depth|inlineStyleRanges|entityRanges|data|2885024|C[i][j]%2B=A[i][k]*B[k][j]|code-block|syntax|javascript|2885025|当k=2时崩溃，这是因为元组A[i]只有2个值，所以你只能在它出错之前调用它到Ai。|offset|length|style|CODE|2885026|编辑:也听听Gerard的答案，你的C是错的。应该是C=[[0+for+row+in+range(len(A))]+for+col+in+range(len(A[0]))]。|2885027|提示:您可以将第一个循环替换为乘法，因此它将是C=[[0]*len(A)+for+col+in+range(len(A[0]))]|2885028|entityMap^0|0|0|E|4|0|Q|1O|0|N|16|0^^$0|@$1|2|3|4|5|6|7|S|8|@]|9|@]|A|$]]|$1|B|3|C|5|D|7|T|8|@]|9|@]|A|$E|F]]|$1|G|3|H|5|6|7|U|8|@$I|V|J|W|K|L]]|9|@]|A|$]]|$1|M|3|N|5|6|7|X|8|@$I|Y|J|Z|K|L]]|9|@]|A|$]]|$1|O|3|P|5|6|7|10|8|@$I|11|J|12|K|L]]|9|@]|A|$]]|$1|Q|3|-4|5|6|7|13|8|@]|9|@]|A|$]]]|R|$]]

The fault occurs here:

<pre><code>C[i][j]+=A[i][k]*B[k][j]
</code></pre>

It crashes when k=2. This is because the tuple <code>A[i]</code> 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 <code>C=[[0 for row in range(len(A))] for col in range(len(A[0]))]</code>.

Just a tip: you could replace the first loop with a multiplication, so it would be <code>C=[[0]*len(A) for col in range(len(A[0]))]</code>

blocks|key|228511|text|这是不正确的初始化。您与col1互换了行！|type|unstyled|depth|inlineStyleRanges|entityRanges|data|228512|C+=+[[0+for+row+in+range(len(A))]+for+col+in+range(len(B[0]))]|code-block|syntax|javascript|228513|正确的初始化应该是|228514|C+=+[[0+for+col+in+range(len(B[0]))]+for+row+in+range(len(A))]|228515|另外，我建议使用更好的命名约定。将在调试方面对您有很大帮助。例如：|228516|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]+%2B=+A[i][k]+*+B[k][j]
++++return+C|228517|你可以做的更多，但是你明白了.|228518|entityMap^0|0|0|0|0|0|0|0^^$0|@$1|2|3|4|5|6|7|S|8|@]|9|@]|A|$]]|$1|B|3|C|5|D|7|T|8|@]|9|@]|A|$E|F]]|$1|G|3|H|5|6|7|U|8|@]|9|@]|A|$]]|$1|I|3|J|5|D|7|V|8|@]|9|@]|A|$E|F]]|$1|K|3|L|5|6|7|W|8|@]|9|@]|A|$]]|$1|M|3|N|5|D|7|X|8|@]|9|@]|A|$E|F]]|$1|O|3|P|5|6|7|Y|8|@]|9|@]|A|$]]|$1|Q|3|-4|5|6|7|Z|8|@]|9|@]|A|$]]]|R|$]]

This is incorrect initialization. You interchanged row with col! 

<pre><code>C = [[0 for row in range(len(A))] for col in range(len(B[0]))]
</code></pre>

Correct initialization would be

<pre><code>C = [[0 for col in range(len(B[0]))] for row in range(len(A))]
</code></pre>

Also I would suggest using better naming conventions. Will help you a lot in debugging. For example:

<pre><code>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
</code></pre>

You can do a lot more, but you get the idea...

blocks|key|229218|text|如果你真的不想使用numpy，你可以这样做：|type|unstyled|depth|inlineStyleRanges|offset|length|style|CODE|entityRanges|data|229219|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)+++++++|code-block|syntax|javascript|229220|结果：|229221|>>>+matmult(x,y)
[[12,+18],+[27,+42],+[42,+66],+[57,+90]]
>>>+mx+*+my
matrix([[12,+18],
++++++++[27,+42],
++++++++[42,+66],
++++++++[57,+90]])|229222|entityMap^0|9|5|0|0|0|0^^$0|@$1|2|3|4|5|6|7|Q|8|@$9|R|A|S|B|C]]|D|@]|E|$]]|$1|F|3|G|5|H|7|T|8|@]|D|@]|E|$I|J]]|$1|K|3|L|5|6|7|U|8|@]|D|@]|E|$]]|$1|M|3|N|5|H|7|V|8|@]|D|@]|E|$I|J]]|$1|O|3|-4|5|6|7|W|8|@]|D|@]|E|$]]]|P|$]]

If you really don't want to use <code>numpy</code> you can do something like this:

<pre><code>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) 
</code></pre>

Result:

<pre><code>&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]])
</code></pre>

blocks|key|2885126|text|下面是我多年前用纯Python编写的矩阵/向量例程的一些简短代码：|type|unstyled|depth|inlineStyleRanges|entityRanges|data|2885127|'''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+%2B+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%2Blhs
++++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+%2B+1.0)
++++++++++++v+/=+t
++++++++return+v,+2.0+*+t**2+/+(sigma+%2B+t**2)
++++def+polyval(+self,+x+):
++++++++'Vec([6,3,4]).polyval(5)+evaluates+to+6*x**2+%2B+3*x+%2B+4+at+x=5'
++++++++return+reduce(+lambda+cum,c:+cum*x%2Bc,+self,+0.0+)
++++def+ratval(+self,+x+):
++++++++'Vec([10,20,30,40,50]).ratfit(5)+evaluates+to+(10*x**2+%2B+20*x+%2B+30)+/+(40*x**2+%2B+50*x+%2B+1)+at+x=5.'
++++++++degree+=+len(self)+/+2
++++++++num,+den+=+self[:degree%2B1],+self[degree%2B1:]+%2B+[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:+%25s+*+%25s'+%25+(`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+%25d+must+match+vector+length+%25d'+%25+(self.rows,+len(b))
++++++++x+=+self._solve(+b+)
++++++++diff+=+b+-+self.mmul(x)
++++++++maxdiff+=+diff.dot(diff)
++++++++for+i+in+range(10):
++++++++++++xnew+=+x+%2B+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%2B1,+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%2B1,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+%25s'+%25+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%2B1)
++++++++++++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%2B1)+*+self[i][i]
++++++++++++++++q,+r+=+(self+-+shift).qr()
++++++++++++++++self+=+r.mmul(q)+%2B+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%2B1:]))+/+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%2B1:]))+):
++++++++++++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%2Bhigh-low)/steps
++++++++xvec+=+[low%2Bh/2.0%2Bh*i+for+i+in+range(steps)]
++++else:
++++++++scale,+base+=+(0.0%2Bhigh-low)/2.0,+(0.0%2Bhigh%2Blow)/2.0
++++++++xvec+=+[base%2Bscale*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+%2B+b*x+%2B+c)/(d*x**2+%2B+e*x+%2B+1)'
++++return+Mat([[x**n+for+n+in+range(degree,-1,-1)]%2B[-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%2Bj%2B1)'
++++return+genmat(m,n,+lambda+i,j:+1.0/(i%2Bj%2B1.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%2B2j,2,3,4])
++++b+=+Table([5,6,7,8])
++++C+=+Table([a,b])
++++print+'a%2Bb',+a%2Bb
++++print+'2%2Ba',+2%2Ba
++++print+'a/5.0',+a/5.0
++++print+'2*a%2B3*b',+2*a%2B3*b
++++print+'a%2BC',+a%2BC
++++print+'3%2BC',+3%2BC
++++print+'C%2Bb',+C%2Bb
++++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+%2B+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)%2B2*math.cos(.7*x%2B.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)%2B2*math.cos(.7*x%2B.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)|code-block|syntax|javascript|2885128|entityMap^0|0|0^^$0|@$1|2|3|4|5|6|7|I|8|@]|9|@]|A|$]]|$1|B|3|C|5|D|7|J|8|@]|9|@]|A|$E|F]]|$1|G|3|-4|5|6|7|K|8|@]|9|@]|A|$]]]|H|$]]

Here's some short and simple code for matrix/vector routines in pure Python that I wrote many years ago:

<pre><code>'''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)
</code></pre>

blocks|key|2885139|text|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]+%2B=+X[i][k]+*+Y[k][j]

for+r+in+result:
+++print(r)|type|code-block|depth|inlineStyleRanges|entityRanges|data|syntax|javascript|2885140|unstyled|entityMap^0|0^^$0|@$1|2|3|4|5|6|7|G|8|@]|9|@]|A|$B|C]]|$1|D|3|-4|5|E|7|H|8|@]|9|@]|A|$]]]|F|$]]

<pre><code>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)
</code></pre>

blocks|key|2885146|text|纯python中的矩阵乘法。|type|unstyled|depth|inlineStyleRanges|entityRanges|data|2885147|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%2B(m1[i][k]*m2[k][j])
++++++++r.append(sums)
++++m.append(r)
++++r=[]
return+m|code-block|syntax|javascript|2885148|entityMap^0|0|0^^$0|@$1|2|3|4|5|6|7|I|8|@]|9|@]|A|$]]|$1|B|3|C|5|D|7|J|8|@]|9|@]|A|$E|F]]|$1|G|3|-4|5|6|7|K|8|@]|9|@]|A|$]]]|H|$]]

Matrix Multiplication in pure python.

<pre><code>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
</code></pre>

blocks|key|229267|text|下面的所有答案都将返回将其转换为矩阵所需的list.Your|type|unstyled|depth|inlineStyleRanges|entityRanges|data|229268|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]+%2B=+A[i][k]+*+B[k][j]

++++C+=+numpy.matrix(C).reshape(len(A),len(B[0]))

++++return+C|code-block|syntax|javascript|229269|entityMap^0|0|0^^$0|@$1|2|3|4|5|6|7|I|8|@]|9|@]|A|$]]|$1|B|3|C|5|D|7|J|8|@]|9|@]|A|$E|F]]|$1|G|3|-4|5|6|7|K|8|@]|9|@]|A|$]]]|H|$]]

All the below answers would return you the list.Your need to convert it to matrix

<pre><code>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
</code></pre>

blocks|key|3099094|text|当我不得不做一些矩阵运算时，我定义了一个新的类来提供帮助。在这样的类中，您可以定义像__add__这样的神奇方法，或者在您的用例中定义__matmul__，允许您定义x+=+a+@+b或a+@=+b而不是matrixMult(a,b)。每个PEP+465在Python3.5中添加了__matmul__。|type|unstyled|depth|inlineStyleRanges|offset|length|style|CODE|entityRanges|data|3099095|我在下面包含了一些实现此功能的代码(我排除了超长的__init__方法，该方法本质上根据传递给它的内容创建一个二维列表self.mat和一个元组self.order+)。|3099096|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]|code-block|syntax|javascript|3099097|注意：|3099098|如果调用了|3099099|__rmatmul__但b不实现__matmul__，则使用|unordered-list-item|3099100|__rmatmul__(例如，如果我想实现2D+list)|3099101|__imatmul__的预乘，则correctly;|3099102|If矩阵在乘法下是不一致的，这意味着它不能相乘，通常是因为它的行数比multiplicand|3099103|中的列数多或少|3099104|3099105|entityMap|0|LINK|mutability|MUTABLE|url|https://www.python.org/dev/peps/pep-0465/^0|16|7|1V|A|2B|9|2L|6|2U|F|3X|A|3C|7|0|0|P|8|1N|8|20|A|0|0|0|0|0|B|C|1|G|A|1|0|B|1|0|B|1|1|0|0^^$0|@$1|2|3|4|5|6|7|1A|8|@$9|1B|A|1C|B|C]|$9|1D|A|1E|B|C]|$9|1F|A|1G|B|C]|$9|1H|A|1I|B|C]|$9|1J|A|1K|B|C]|$9|1L|A|1M|B|C]]|D|@$9|1N|A|1O|1|1P]]|E|$]]|$1|F|3|G|5|6|7|1Q|8|@$9|1R|A|1S|B|C]|$9|1T|A|1U|B|C]|$9|1V|A|1W|B|C]]|D|@]|E|$]]|$1|H|3|I|5|J|7|1X|8|@]|D|@]|E|$K|L]]|$1|M|3|N|5|6|7|1Y|8|@]|D|@]|E|$]]|$1|O|3|P|5|6|7|1Z|8|@]|D|@]|E|$]]|$1|Q|3|R|5|S|7|20|8|@$9|21|A|22|B|C]|$9|23|A|24|B|C]|$9|25|A|26|B|C]]|D|@]|E|$]]|$1|T|3|U|5|S|7|27|8|@$9|28|A|29|B|C]]|D|@]|E|$]]|$1|V|3|W|5|S|7|2A|8|@$9|2B|A|2C|B|C]]|D|@]|E|$]]|$1|X|3|Y|5|S|7|2D|8|@]|D|@]|E|$]]|$1|Z|3|10|5|6|7|2E|8|@]|D|@]|E|$]]|$1|11|3|-4|5|6|7|2F|8|@]|D|@]|E|$]]|$1|12|3|-4|5|6|7|2G|8|@]|D|@]|E|$]]]|13|$14|$5|15|16|17|E|$18|19]]]]

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 <code>__add__</code>, or, in your use-case, <code>__matmul__</code>, allowing you to define <code>x = a @ b</code> or <code>a @= b</code> rather than <code>matrixMult(a,b)</code>. <code>__matmul__</code> was added in Python 3.5 per <a href="https://www.python.org/dev/peps/pep-0465/" rel="nofollow noreferrer">PEP 465</a>.

I have included some code which implements this below (I excluded the prohibitively long <code>__init__</code> method, which essentially creates a two-dimensional list <code>self.mat</code> and a tuple <code>self.order</code> according to what is passed to it)

<pre><code>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]
</code></pre>

Note:

<ul>
<li><code>__rmatmul__</code> is used if <code>b @ a</code> is called and <code>b</code> does not implement <code>__matmul__</code> (e.g. if I wanted to implement premultiplying by a 2D list)</li>
<li><code>__imatmul__</code> is required for <code>a @= b</code> to work correctly;</li>
<li>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</li>
</ul>

blocks|key|2885217|text|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]+%2B=+A[i][k]*B[k][j]
++++return+C|type|code-block|depth|inlineStyleRanges|entityRanges|data|syntax|javascript|2885218|在第二行，您应该更改|unstyled|2885219|C+=+[[0+for+row+in+range(len(B[0]))]+for+col+in+range(len(A))]|2885220|entityMap^0|0|0|0^^$0|@$1|2|3|4|5|6|7|K|8|@]|9|@]|A|$B|C]]|$1|D|3|E|5|F|7|L|8|@]|9|@]|A|$]]|$1|G|3|H|5|6|7|M|8|@]|9|@]|A|$B|C]]|$1|I|3|-4|5|F|7|N|8|@]|9|@]|A|$]]]|J|$]]

<pre><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
</code></pre>

at second line you should change

<pre><code>C = [[0 for row in range(len(B[0]))] for col in range(len(A))]
</code></pre>

blocks|key|230883|text|一个线条：|type|unstyled|depth|inlineStyleRanges|entityRanges|data|230884|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
++++]|code-block|syntax|javascript|230885|解释：|230886|zip(*m2)+-从第二个矩阵中获取一列|230887|zip(m1_r，m2_c)+-从m1行和m2列创建元组|230888|sum(...)+-sum乘法行*列|230889|测试：|230890|m1+=+[[1,+2,+3],+[4,+5,+6]]
m2+=+[[7,+8],+[9,+10],+[11,+12]]
result+=+dot(m1,+m2)
assert+result+==+[[58,+64],+[139,+154]]|230891|entityMap^0|0|0|0|0|0|0|0|0^^$0|@$1|2|3|4|5|6|7|U|8|@]|9|@]|A|$]]|$1|B|3|C|5|D|7|V|8|@]|9|@]|A|$E|F]]|$1|G|3|H|5|6|7|W|8|@]|9|@]|A|$]]|$1|I|3|J|5|6|7|X|8|@]|9|@]|A|$]]|$1|K|3|L|5|6|7|Y|8|@]|9|@]|A|$]]|$1|M|3|N|5|6|7|Z|8|@]|9|@]|A|$]]|$1|O|3|P|5|6|7|10|8|@]|9|@]|A|$]]|$1|Q|3|R|5|D|7|11|8|@]|9|@]|A|$E|F]]|$1|S|3|-4|5|6|7|12|8|@]|9|@]|A|$]]]|T|$]]

One liner:

<pre><code>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
 ]
</code></pre>

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:

<pre><code>m1 = [[1, 2, 3], [4, 5, 6]]
m2 = [[7, 8], [9, 10], [11, 12]]
result = dot(m1, m2)
assert result == [[58, 64], [139, 154]]
</code></pre>

I'm trying to multiply two matrices together using pure Python. Input (<code>X1</code> is a 3x3 and <code>Xt</code> is a 3x2):

<pre><code>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)]
</code></pre>

where Xt is the zip transpose of another matrix. Now here is the code:

<pre><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
</code></pre>

The error that python gives me is this: 

<blockquote>
 IndexError: list index out of range. 
</blockquote>

Now I'm not sure if <code>Xt</code> is recognised as an matrix and is still a list object, but technically this should work.

Matrix Multiplication in pure Python?

Python

我正在尝试使用纯Python将两个矩阵相乘。输入(X1为3x3，Xt为3x2)：X1 =  [[1.0016, 0.0, -16.0514],        [0.0, 10000.0, -40000.0],        [-16.0514, -40000.0, 160513.6437]]Xt =  [(1.0, 1...

问纯Python中的矩阵乘法？
EN

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

问纯Python中的矩阵乘法？EN

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

问纯Python中的矩阵乘法？
EN