Opened 8 years ago

Closed 5 years ago

#10815 closed enhancement (wontfix)

Matrix multiplication and power are very slow

Reported by: frieux Owned by: jason, was
Priority: minor Milestone: sage-duplicate/invalid/wontfix
Component: linear algebra Keywords: matrix, multiplication, slow
Cc: tmonteil Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #15944 Stopgaps:

Description

It seems that the multiplication and the power of a matrix is very slow:

size = 500
A = matrix(RR, size, size)
for i in range(size):
    A[i,i] = random()
    for j in range(8):
        A[i,randint(0,size-1)] = random()

timeit('A*A')

gives

CPU times: user 60.72 s, sys: 0.11 s, total: 60.83 s
Wall time: 61.10 s
500 x 500 dense matrix over Real Field with 53 bits of precision
sage: timeit('A*A')
5 loops, best of 3: 60.6 s per loop

Trying with the 100th power of a 2000*2000 matrix is impossible, according to a colleague working for linbox, taking the 100th power of such a matrix would take no more that 10 seconds.

Is the interface with atlas working well ?

Change History (9)

comment:1 Changed 8 years ago by dsm

I think the issue is that you're doing the computations in RR, the Real Field with 53 bits, rather than in RDF, the Real Double Field.

for size in [5, 10, 20, 50, 100, 200, 500, 1000, 2000]:
    A = matrix(RDF, size, size, dense=True)
    for i in range(size):
        A[i,i] = random()
        for j in range(8):
            A[i,randint(0,size-1)] = random()
    print 'size:', size
    print 'time for one multiplication:'
    timeit('A*A')
    print 'time for A**100:'
    timeit('A**100')

gives (on my notebook)

size: 5
time for one multiplication:
625 loops, best of 3: 79 µs per loop
time for A**100:
625 loops, best of 3: 646 µs per loop
size: 10
time for one multiplication:
625 loops, best of 3: 79.3 µs per loop
time for A**100:
625 loops, best of 3: 648 µs per loop
size: 20
time for one multiplication:
625 loops, best of 3: 83.3 µs per loop
time for A**100:
625 loops, best of 3: 687 µs per loop
size: 50
time for one multiplication:
625 loops, best of 3: 154 µs per loop
time for A**100:
625 loops, best of 3: 1.26 ms per loop
size: 100
time for one multiplication:
625 loops, best of 3: 499 µs per loop
time for A**100:
125 loops, best of 3: 4.11 ms per loop
size: 200
time for one multiplication:
125 loops, best of 3: 3.24 ms per loop
time for A**100:
25 loops, best of 3: 26.2 ms per loop
size: 500
time for one multiplication:
25 loops, best of 3: 35 ms per loop
time for A**100:
5 loops, best of 3: 281 ms per loop
size: 1000
time for one multiplication:
5 loops, best of 3: 237 ms per loop
time for A**100:
5 loops, best of 3: 1.9 s per loop
size: 2000
time for one multiplication:
5 loops, best of 3: 1.64 s per loop
time for A**100:
5 loops, best of 3: 13.4 s per loop

which I think is pretty close to your colleague's expectation.

comment:2 Changed 6 years ago by kcrisman

  • Status changed from new to needs_info
  • Type changed from defect to enhancement

So... what should we do here? I think that maybe the best thing is to find more ways in various reference and tutorial sources to say "use RDF for speed", but I don't know how realistic that is. Anyway, not a bug, so setting to "needs info".

comment:3 Changed 6 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:4 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:5 Changed 5 years ago by mmezzarobba

  • Milestone changed from sage-6.2 to sage-duplicate/invalid/wontfix
  • Priority changed from critical to minor
  • Status changed from needs_info to needs_review

comment:6 Changed 5 years ago by rws

May also be effected by #15912.

comment:7 Changed 5 years ago by tmonteil

  • Dependencies set to #15944

@rws: this is not related to printing, but to the fact that RR privides emulated floating-point arithmetics (by MPFR), which is inherently slower than CPU implementation provided by RDF. Moreover, some optimized libraries (like atlas) require CPU doubles to work with, and are not used in such a case.

I was the one surprised by Frederic's timings. As @kcrisman, i think there is still an issue about misleadingly using RR as the default implementation of real floating point (which seems to come from its generic name).

This documentation issue should be fixed by the forthcomig tutorial #15944.

comment:8 Changed 5 years ago by rws

  • Status changed from needs_review to positive_review

comment:9 Changed 5 years ago by vbraun

  • Resolution set to wontfix
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.