On 2014-12-07, Jernej <azi.std...@gmail.com> wrote:
> ------=_Part_3308_1073346028.1417954673282
> Content-Type: multipart/alternative;
> boundary="----=_Part_3309_733817265.1417954673282"
>
> ------=_Part_3309_733817265.1417954673282
> Content-Type: text/plain; charset=UTF-8
>
> Just for fun I wrote an equivalent program in C and tested the Sage
> function and an equivalent C program on 1 instance of the problem.
>
>===
> $ time ./a.out test
>
> real 0m0.366s
> user 0m0.306s
> sys 0m0.060s
>===
>
> And Sage
>
>===
> $ time sage test.sage test
>
> real 2m6.285s
> user 2m5.256s
> sys 0m0.858s
>===
>
> I thought that may be interesting to you as well.
>
> Best,
>
> Jernej
>
>
> On Wednesday, 3 December 2014 22:33:46 UTC+1, Jernej wrote:
>>
>> Dear sage-support,
>>
>> I have stumbled into a performance bottleneck in one of my Sage programs.
>> I would like to share the relevant problem here in hope anyone has a
>> constructive suggestion for optimizing the given program.
>>
>> I am given a n x n, (0,1) matrix C where n < 20. C has up to 30% of
>> nonzero elements.
>>
>> I define the inner product <u,v> = u (2I - C) v^t and need to compute the
>> set S of all (0,1) n-vectors u such that <u,u> = 2 and <u,j> = -1, where j
>> is the vector with all ones. Finally I need to record all pairs u,v from S
>> such that <u,v> == -1 or 0.
>>
>> The one optimization that I use is to cache the product u*(2I-C). Other
>> than that I don't see any other way to optimize the program hence it looks
>> like
<u,j>=-1 is a linear condition on u.
Shouldn't you only pick up u from the corresponding hyperplane?
>>
>> ====
>> global vectors
>> field = RR # this looks like the fastest option
RDF should be faster, IMHO.
>>
>> def init_vectors(n):
>> global vectors
>> vectors = [matrix(field,b) for b in product(*[[0,1]]*n)]
>> def foo(C,output):
>>
>> global vectors
>>
>> n = C.nrows()
>> # we compute the inverse in QQ just to gain a bit of numerical
>> stability
>> D = (2*identity_matrix(n)-C).inverse()
>> D = Matrix(field,D)
>>
>> j = matrix(field, tuple([1 for _ in xrange(n)]))
>> cur = 1
>>
>> vec2int = {}
>> cache = {}
>>
>> for bm in vectors:
>> val = bm*D
>> if abs( (val*bm.transpose())[0,0]-2) <= 0.000001 and
>> abs((val*j.transpose())[0,0]+1) <= 0.000001:
>> vec2int[cur] = bm
>> output.write(str(el for el in bm[0])+'\n')
>> cache[cur] = val
>> cur+=1
>>
>> for i in xrange(1, cur):
>> for j in xrange(i+1, cur):
>> iv = (cache[i]*vec2int[j].transpose())[0,0]
>> if abs(iv) <= 0.0000001 or abs(iv+1) <= 0.000001:
>> output.write(str(i) + ' ' + str(j) + '\n')
>> ===
>>
>> For convenience's sake I have attached the result of %prun on one instance
>> of the problem. From the output of prun I suppose it would make sense to
>> cache the transpose of the vectors as well, though I am not sure this is
>> going to make the crucial difference.
>>
>> On average it takes 200s to process one matrix and given that I have
>> millions of such matrices it would take years to process all the input.
>>
>> Other than rewriting the whole thing in C I currently do not see any
>> viable solution. Hence I am wondering: Do you guys happen to see any clever
>> optimization? Is there a way to compute the named product more efficiently?
You could use Cython instead of C, and numpy arrays instead of Sage
matrices (numpy arrays interface quite well with Cython).
Alas, currently plain Sage matrices are very slow.
As your matrices are quite small, alternatively you might experiment
with computing modulo p, for p a prime or order bigger than n^2.
>>
>> Best,
>>
>> Jernej
>>
>
--
You received this message because you are subscribed to the Google Groups
"sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sage-support+unsubscr...@googlegroups.com.
To post to this group, send email to sage-support@googlegroups.com.
Visit this group at http://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.
