Repeated dot products
Hi, Suppose I have an array of shape: (n, k, k). In this case, I have n k-by-k matrices. My goal is to compute the product of a (potentially large) user-specified selection (with replacement) of these matrices. For example, x = [0,1,2,1,3,3,2,1,3,2,1,5,3,2,3,5,2,5,3,2,1,3,5,6] says that I want to take the 0th matrix and dot it with the 1st matrix and dot that product with the 2nd matrix and dot that product with the 1st matrix again and so on... Essentially, I am looking for efficient ways of doing this. It seems like what I *need* is for dot to be a ufunc with a reduce() operator. Then I would construct an array of the matrices, as specified by the input x. For now, I am using a python loop and this is unbearable.
prod = np.eye(k) for i in x: ... prod = dot(prod, matrices[i]) ...
Is there a better way?
On Sat, Dec 12, 2009 at 4:55 PM, T J <tjhnson@gmail.com> wrote:
Hi,
Suppose I have an array of shape: (n, k, k). In this case, I have n k-by-k matrices. My goal is to compute the product of a (potentially large) user-specified selection (with replacement) of these matrices. For example,
x = [0,1,2,1,3,3,2,1,3,2,1,5,3,2,3,5,2,5,3,2,1,3,5,6]
says that I want to take the 0th matrix and dot it with the 1st matrix and dot that product with the 2nd matrix and dot that product with the 1st matrix again and so on...
Essentially, I am looking for efficient ways of doing this. It seems like what I *need* is for dot to be a ufunc with a reduce() operator. Then I would construct an array of the matrices, as specified by the input x. For now, I am using a python loop and this is unbearable.
prod = np.eye(k) for i in x: ... prod = dot(prod, matrices[i]) ...
Is there a better way?
I don't know about numpy, but I was using python reduce for similar:
d = np.eye(2) li = [d,d,d,d,d,d,d,d] reduce(np.dot,li) array([[ 1., 0.], [ 0., 1.]])
lii = [1,2,1,2,3] reduce(lambda x,i : np.dot(x,li[i]),lii) array([[ 1., 0.], [ 0., 1.]])
quickly written from memory, not sure its always correct, and no idea about speed compared to loop Josef
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Le samedi 12 décembre 2009 à 01:55 -0800, T J a écrit :
Is there a better way?
You may have a look at http://scipy.org/Cookbook/MultiDot Several alternatives are proposed. Cheers, -- Denis
2009/12/12 T J <tjhnson@gmail.com>:
Hi,
Suppose I have an array of shape: (n, k, k). In this case, I have n k-by-k matrices. My goal is to compute the product of a (potentially large) user-specified selection (with replacement) of these matrices. For example,
x = [0,1,2,1,3,3,2,1,3,2,1,5,3,2,3,5,2,5,3,2,1,3,5,6]
says that I want to take the 0th matrix and dot it with the 1st matrix and dot that product with the 2nd matrix and dot that product with the 1st matrix again and so on...
Essentially, I am looking for efficient ways of doing this. It seems like what I *need* is for dot to be a ufunc with a reduce() operator. Then I would construct an array of the matrices, as specified by the input x. For now, I am using a python loop and this is unbearable.
prod = np.eye(k) for i in x: ... prod = dot(prod, matrices[i]) ...
Is there a better way?
You are right that numpy/scipy should have a matrix product ufunc. In fact it should have ufunc versions of all the linear algebra tools. It even has the infrastructure (generalized ufuncs) to support such a thing in a very general way. Sadly this infrastructure is basically unused. Your best workaround depends on the sizes of the various objects involved. If k is large enough, then a k by k by k matrix multiply will be slow enough that it doesn't really matter what else you do. If x is much longer than n, you will want to avoid constructing a len(x) by k by k array. If n is really small and x really large, you might even win by some massively clever Knuthian operation that found repeated substrings in x and evaluated each corresponding matrix product only once. If on the other hand len(x) is much shorter than n, you could perhaps benefit from forming an intermediate len(x) by k by k array. This array could then be used in a vectorized reduce operation, if we had one. Normally, one can work around the lack of a vectorized matrix multiplication by forming an (n*k) by k matrix and multiplying it, but I don't think this will help any with reduce. If k is small, you can multiply two k by k matrices by producing the k by k by k elementwise product and then adding along the middle axis, but I don't think this will work well with reduction either. So I have just two practical solutions, really: (a) use python reduce() and the matrix product, possibly taking advantage of the lack of need to produce an intermediate len(x) by k by k matrix, and live with python in that part of the loop, or (b) use the fact that in recent versions of numpy there's a quasi-undocumented ufuncized matrix multiply built as a test of the generalized ufunc mechanism. I have no idea whether it supports reduce() (and if it doesn't adding it may be an unpleasant experience). Good luck, Anne P.S. if you come up with a good implementation of the repeated substring approach I'd love to hear it! -A
On 12/12/2009 22:55, T J wrote:
Hi,
Suppose I have an array of shape: (n, k, k). In this case, I have n k-by-k matrices. My goal is to compute the product of a (potentially large) user-specified selection (with replacement) of these matrices. For example,
x = [0,1,2,1,3,3,2,1,3,2,1,5,3,2,3,5,2,5,3,2,1,3,5,6]
TJ, what are your n, k, len(x) ? _dotblas.dot is fast: dot( 10x10 matrices ) takes ~ 22 usec on my g4 ppc, which is ~ 15 clock cycles (700 MHz) per mem access * +. A hack to find repeated pairs (or triples ...) follows. Your sequence above has only (3,2) 4 times, no win. (Can someone give a probabilistic estimate of the number of non-overlapping pairs in N letters from an alphabet of size A ?) #!/usr/bin/env python # numpy-discuss 2009 12dec TJ repeated dot products from __future__ import division from collections import defaultdict import numpy as np __version__ = "2010 7jan denis" def pairs( s, Len=2 ): """ repeated non-overlapping pairs (substrings, subwords) "abracadabra" -> ab ra [[0 7] [2 9]], not br Len=3: triples, 4 ... """ # bruteforce # grow repeated 2 3 ... ? pairs = defaultdict(list) for j in range(len(s)-Len+1): pairs[ s[j:j+Len] ].append(j) min2 = filter( lambda x: len(x) > 1, pairs.values() ) min2.sort( key = lambda x: len(x), reverse=True ) # remove overlaps -- # (if many, during init scan would be faster) runs = np.zeros( len(s), np.uint8 ) run = np.ones( Len, np.uint8 ) run[0] = Len chains = [] for ovchain in min2: chain = [] for c in ovchain: if not runs[c:c+Len].any(): runs[c:c+Len] = run chain.append(c) if len(chain) > 1: chains.append(chain) return (chains, runs) #............................................................................... if __name__ == "__main__": import sys abra = "abracadabra" alph = 5 randlen = 100 randseed = 1 exec( "\n".join( sys.argv[1:] )) # Test= ... print "pairs( %s ) --" % abra print pairs( abra ) # ab [0, 7], br [2, 9]] print pairs( abra, 3 ) # abr [0, 7] np.random.seed( randseed ) r = np.random.random_integers( 1, alph, randlen ) chains, runs = pairs( tuple(r) ) npair = sum([ len(c) for c in chains ]) print "%d repeated pairs in %d random %d" % (npair, randlen, alph) # 35 repeated pairs in 100 random 5 (prob estimate this ?) # 25 repeated pairs in 100 random 10
participants (5)
-
Anne Archibald -
denis -
Denis Laxalde -
josef.pktd＠gmail.com -
T J