matAtBA computes the matrix expression $$A^tBA$$.

matABAt computes the matrix expression $$ABA^t$$.

colmat2vec returns the vector for 1 column matrix.

rowmat2vec returns the vector for a 1 rowed matrix.

Any requests for features should be addressed to the author, [ernesto . adorio @ gmail . com], remove the internal spaces.

# -*- coding: utf-8 -*- """ file matlib.py author Ernesto P. Adorio UPDEPP (UP Clark), Pampanga ernesto.adorio@gmail.com revisions Ver 0.0.1 initial release 2009.01.16 added matdots,matrandom, isiterable Ver 0.0.2 2009.10.12 added vec2colmat,vec2rowmat, mat2vec, matSelectCols, matDelCols, matInsertConstCol,matreadstring Ver 0.0.3 2009.12.22 added matSubmat, matRows, matCols, matmatt, matxtx, matvectmatvec,matsplitbycolumn revised matvec Ver 0.0.4 2010.06.20 added matAtBA, matABAt.colmat2vec, rowmat2vec """ def matreadstring(s, sep = " "): """ Each line must be complete row. comment lines start with a #. version 0.0.1 no varnames yet! """ M = [] for line in s.split("\n"): if line: if not line.startswith("#"): M.append([float(x) for x in (line.split(sep))]) return M def matsplitbycolumn(X, col=-1): """ Splits a matrix X by column. """ if col == -1 or col == len(X[0])-1: # X2 will be a vector not a matrix. X1 = [x[:-1] for x in X] X2 = [x[-1] for x in X] else: #X2 is a matrix. X1 = [x[:col] for x in X] X2 = [X[col:-1] for x in X] return X1, X2 def tabular(hformat, headers, bformat, M): # added nov.29, 2008. # prints the table data. nrows = len(M) ncols = len(M[0]) # print headers. for j, heading in enumerate(headers): print hformat[j] % heading, print # print the body. for i, row in enumerate(M): for j, col in enumerate(M[i]): print bformat[j] % M[i][j], print print def vecadd(X,Y): n = len(X) if n != len(Y): return None return [x + y for x,y in zip(X,Y)] def vecsub(X, Y): n = len(X) if n != len(Y): raise ArgumentError, "incompatible vectors in vecsub." return [x - y for x,y in zip(X,Y)] def eye(m, n= None): if n is None: n = m B= [[0]* n for i in range(m)] for i in range(m): B[i][i] = 1.0 return B matIdentity = eye matiden = eye iden = eye def vec2colmat(X): """ Retuns a 1 column matrix out of array or vector X. """ return ([[x] for x in X]) def vec2rowmat(X): """ Retuns a 1 row matrix out of array or vector X. """ return ([[x for x in X]]) def mat2vec(M, column=0): """ Returns a vector from a column of matrix M. """ try: M[0][column] return [m[column] for m in M] except: return M def colmat2vec(M): """ Transform a column matrix to a row vector. """ if len(M[0]) == 1: return [x[0] for x in M] def rowmat2vec(M): """ Transform a row matrix to a row vector. """ if len(M) == 1: return([x for x in M[0]]) def matvectmatvec(x,M): """ Computes x^t M x where x is a column vector. Result is a scalar. """ return dot(x,matvec(M,x)) def matSelectCols(M, jindices): """ Extracts a submatrix from M with col indices in jindices. Negative indices are properly handled too by Python, no need for adjustment. """ N = [] n = len(M) for i in range(n): N.append([M[i][j] for j in jindices]) return N matCols = matSelectCols def matRows(M, iindices): """ Returns selected rews of M as a matrix. """ return [M[i:] for i in iindices] def matDelCols(M, colindices): """ M - input matrix colindices = array indices of columns to delete. Returns matrix from M with columns deleted. """ ncols = len(M[0]) nindices = len(colindices) # Adjust for negative indices. for j in range(nindices): if colindices[j] < 0: colindices[j] += ncols jindices = [] for j in range(ncols): if j not in colindices: jindices.append(j) return matSelectCols(M, jindices) def matSubmat(M, rindices, jindices): # Returns a submatrix selected from M. N = [] for i in rindices: N.append([M[i][j] for j in jindices]) return N submatrix = matSubmat def matInsertConstCol(X, column, c = 1.0, inplace= True): """ Inserts a constant column to vector or matrix X at position column. NEVER forget that indicing starts at zero. """ if not inplace: Xcopy = [x[:] for x in X] else: Xcopy = X for i in range(len(Xcopy)): if type(Xcopy[i]) == type(0.0): Xcopy[i] = [c, X[i]] else: Xcopy[i].insert(column, c) return Xcopy def matxtx(X): # returns the matrix of the coefficients of the # normal equations for least squares computation. # This should be more efficient than calling a # matrix multiplication on t(X) and X. # same as function call mattmat(X, X) m, n = matdim(X) if m and n== 0: return sum([x for x in X]) M = [[0.0]* n for i in range(n)] for i in range(n): for j in range(i, n): dot = 0.0 for r in range(m): dot += X[r][i] * X[r][j] M[i][j] = dot if i != j: M[j][i] = dot return M xtx = matxtx def matzero(m, n = None): """ Returns an m by n zero matrix. """ if n is None: n = m return [[0]* n for i in range(m)] def matdiag(D): """ Returns a diagonal matrix with diagonal elements from D. """ n = len(D) A = [[0] * n for i in range(n)] for i in range(n): A[i][i] = D[i] return A def diag(A): """ Returns diagonal elements of A as a vector """ #return [A[i][i] for i in range(len(A))] # Here is a version which returns a vector. return [A[i][i] for i in range(len(A))] def matcol(X, j): # Returns the jth column of matrix X. nrows = len(X) return [X[i][j] for i in range(nrows)] def trace(A): """ Returns the trace of a matrix. """ return sum([A[i][i] for i in range(len(A))]) def matadd(A, B): """ Returns C = A + B. """ try: m = len(A) if m != len(B): return None n = len(A[0]) if n != len(B[0]): return None C = matzero(m, n) for i in range(m): for j in range(n): C[i][j] = A[i][j] + B[i][j] return C except: return None def matsub(A, B): """ returns C = A - B. """ try: m = len(A) if m != len(B): return None n = len(A[0]) if n != len(B[0]): return None C = matzero(m, n) for i in range(m): for j in range(n): C[i][j] = A[i][j] - B[i][j] return C except: return None def matcopy(A): #print "Inside matcopy:", A B = [] for a in A: #print a try: B.append(a[:]) except: B.append(a) return B def matkmul(A, k): """ Multiplies each element of A by k. """ B = matcopy(A) for i in range(len(A)): for j in range(len(A[0])): B[i][j] *= k return B def transpose(A): """ Returns the transpose of A. """ m,n = matdim(A) At = [[0] * m for j in range(n)] for i in range(m): for j in range(n): At[j][i] = A[i][j] m,n = matdim(At) return At matt = transpose mattrans = transpose tr = transpose def matdim(A): # Returns the number of rows and columns of A. m = n = 0 try: iter(A) m = len(A) except: pass try: iter(A[0]) n = len(A[0]) except: pass return (m, n) def matprod(A, B,debug= False): """ Computes the product of two matrices. 2009.01.16 Revised for matrix or vector B. A and B are matrices. If one of them is a vector, it must be transformed into a matrix with one row or one column. """ m, n = matdim(A) p, q = matdim(B) if n!= p: ErrMsg = "matprod: Incompatible matrices, n1 , m2= %s,%s" % (len(A[0]) , len(B)) raise ValueError, ErrMsg try: if iter(B[0]): q = len(B[0]) except: q = 1 if debug: print "inside matprod(),A,B matrices::" matprint(A) matprint(B) print "dimensions:", m,n, p,q print "inside matprod(), multiplying:" C = matzero(m, q) for i in range(m): for j in range(q): t = sum([A[i][k] * B[k][j] for k in range(p)]) C[i][j] = t return C matmul = matprod def matvec(A, y): """ Returns the product of matrix A with vector y. Revision: dec. 22, 2009: this version should work with one column matrices y. """ m = len(A) n = len(A[0]) try: y[0][0] out = [0] * m for i in range(m): for j in range(n): out[i] += A[i][j] * y[j][0] return out except: out = [0] * m for i in range(m): for j in range(n): out[i] += A[i][j] * y[j] return out def mattvec(A, y): """ Returns the vector A^t y. """ At = transpose(A) return matvec(At, y) def dot(X, Y): """ Dot product of vectors X and Y. """ return sum(x* y for (x,y) in zip(X,Y)) def matdots(X): # Added Jan 16, 2009. # Returns the matrix of dot products of the column vectors # This is the same as X^t X. (nrow, ncol) = matdim(X) M = [[0.0] * ncol for i in range(ncol)] for i in range(ncol): for j in range(i+1): dot = sum([X[p][i]* X[p][j] for p in range(ncol)]) M[i][j] = dot if i != j: M[j][i] = M[i][j] return M def matABAt(A,B): """ Returns the product ABA^t. """ AB = matmul(A,B) return matmul(AB, transpose(A)) def matAtBA(A,B): """ Returns the product A^t B A. """ AtB = matmul(transpose(A), B) return matmul(AtB, A) def mattmat(A, B): """ Returns the product [transpose(A) B] if B = A, use matxtx instead. """ AtB = matprod(transpose(A), B) return AtB def matmatt(A,B): """ Returns A B^t added dec,22,2009 """ return matmul (A, mattrans(B)) def matrandom(nrow, ncol = None): # Added Jan. 16, 2009 if ncol is None: ncol = nrow R = [] for i in range(nrow): R.append([random.random() for j in range(ncol)]) return R def matunitize(X, inplace = False): # Added jan. 16, 2009 # Transforms each vector in X to have unit length. if not inplace: V = [x[:] for x in X] else: V = X nrow = len(X) ncol = len(X[0]) for j in range(ncol): recipnorm = sum([X[j][j]**2 for j in range(ncol)]) for i in range(nrow): V[i][j] *= recipnorm return V def matprint(A,format= "%8.3f"): #prints the matrix A using format m, n = matdim(A) if m: for i,row in enumerate(A): if n: for c in row: print format % c, print else: print row else: if n==0: print A else: for c in A: print format %A print def mataugprint(A,Y, format= "%8.4f"): #prints the augmented matrix A|Y using format try: ycols = len(Y[0]) except: ycols = 1 for i,row in enumerate(A): for c in row: print format % c, print "|", if ycols == 1: print format % Y[i] else: for y in Y[i]: print format % Y[i], print def gjinv(AA,inplace = False): """ Determines the inverse of a square matrix BB by Gauss-Jordan reduction. """ n = len(AA) B = eye(n) if not inplace: A = [row[:] for row in AA] else: A = AA for i in range(n): #Divide the ith row by A[i][i] m = 1.0 / A[i][i] for j in range(i, n): A[i][j] *= m # # this is the same as dividing by A[i][i] for j in range(n): B[i][j] *= m #lower triangular elements. for k in range(i+1, n): m = A[k][i] for j in range(i+1, n): A[k][j] -= m * A[i][j] for j in range(n): B[k][j] -= m * B[i][j] #upper triangular elements. for k in range(0, i): m = A[k][i] for j in range(i+1, n): A[k][j] -= m * A[i][j] for j in range(n): B[k][j] -= m * B[i][j] return B matinverse = gjinv matinv = gjinv inverse = gjinv inv = gjinv def Test(): X = [1,1,1] print dot(X, X) AA = [[1,2,3], [4,5,8], [9,7,6]] BB = eye(3) print "Identity matrix eye(3):" matprint(BB) print "inputs:" print AA print "product" matprint(matprod(AA, AA)) print "inverse of AA:" BB = gjinv(AA) matprint(BB) print "product of AA and its inverse:" matprint(matprod(AA ,BB)) if __name__ == "__main__": Test()

## No comments:

## Post a Comment