Opened 9 years ago

Closed 6 years ago

# Matrix multiplication and power are very slow

Reported by: Owned by: frieux jason, was minor sage-duplicate/invalid/wontfix linear algebra matrix, multiplication, slow tmonteil N/A #15944

### 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 ?

### comment:1 Changed 9 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 7 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 7 years ago by jdemeyer

• Milestone changed from sage-5.11 to sage-5.12

### comment:4 Changed 6 years ago by vbraun_spam

• Milestone changed from sage-6.1 to sage-6.2

### comment:5 Changed 6 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 6 years ago by rws

May also be effected by #15912.

### comment:7 Changed 6 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 6 years ago by rws

• Status changed from needs_review to positive_review

### comment:9 Changed 6 years ago by vbraun

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