Matrix_generic_dense.det() is quite slow compared to Numpy

[陳孜穎]

I am doing a numerical analysis which needs to calculate many determinants of "huge" dense matrices. The sizes of the matrix are less than 100*100.

I found the computation of the determinant is quite slow with Matrix_generic_dense. I followed the method on What is the time complexity of numpy.linalg.det? and made minor modifications to the script to benchmark for Sagemath.
1. np.arange(1,10001, 100) => np.arange(1, 292, 10) because the speed is too slow, I have to reduce the testing range to make the script finish within a reasonable time.
2. np.ones((size, size))  => ones_matrix(RR, size, size)
3. timeit('np.linalg.det(A)', globals={'np':np, 'A':A}, number=1) => timeit('A.det(algorithm="df")', globals={'A': A}, number=1)

Here is the full script I used.

from timeit import timeit

import matplotlib.pyplot as plt
import numpy as np
from sage.rings.real_mpfr import RR
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

sizes = np.arange(1, 292, 10)
times = []

for size in sizes:
    A = ones_matrix(RR, size, size)
    time = timeit('A.det(algorithm="df")', globals={'A': A}, number=1)
    times.append(time)
    print(size, time)

sizes = sizes.reshape(-1, 1)
times = np.array(times).reshape(-1, 1)

quad_sizes = PolynomialFeatures(degree=2).fit_transform(sizes)
quad_times = LinearRegression().fit(quad_sizes, times).predict(quad_sizes)
cubic_sizes = PolynomialFeatures(degree=3).fit_transform(sizes)
cubic_times = LinearRegression().fit(cubic_sizes, times).predict(cubic_sizes)
quartic_sizes = PolynomialFeatures(degree=4).fit_transform(sizes)
quartic_times = LinearRegression().fit(quartic_sizes,
                                       times).predict(quartic_sizes)

plt.scatter(sizes, times, label='Data', color='k', alpha=0.5)
plt.plot(sizes, quad_times, label='Quadratic', color='r')
plt.plot(sizes, cubic_times, label='Cubic', color='g')
plt.plot(sizes, quartic_times, label='Quartic', color='b')
plt.xlabel('Matrix Dimension n')
plt.ylabel('Time (seconds)')
plt.legend()
plt.show()

This is the time complexity diagram on the linked post. The data from the testing with Numpy. It only took about 20 seconds to calculate a 10000*10000 matrix.

This is the testing with matrix(RR) with the df algorithm.

I also tested with the "hessenberg" algorithm which got worse.

Compared to Numpy, Sagemath spent 800 seconds calculating a 300*300 matrix. It seems that Sagemath still runs the algorithm under O(n^3) (the curves of the cubic model and the quartic model overlap) but with a large constant. Why Sagemath is so slow? Is it possible to speed up?

[John H Palmieri]

One way to speed it up might be to work with RDF or QQ or RLF. The documentation for the generic determinant method says, "Note that for matrices over most rings, more sophisticated algorithms can be used."

sage: %time ones_matrix(RDF, 600, 600).determinant()
CPU times: user 78.6 ms, sys: 7.6 ms, total: 86.2 ms
Wall time: 77.9 ms
0.0
sage: %time ones_matrix(QQ, 600, 600).determinant()
CPU times: user 280 ms, sys: 36.3 ms, total: 317 ms
Wall time: 325 ms
0
sage: %time ones_matrix(RLF, 600, 600).determinant()
CPU times: user 32.7 s, sys: 183 ms, total: 32.9 s
Wall time: 40.3 s
0

[陳孜穎]

Is there any space to improve the computation of MPFR which is used by RealField?

John H Palmieri 在 2023年1月6日 星期五凌晨3:01:37 [UTC+8] 的信中寫道：

[John H Palmieri]

I am not at all an expert in numerical linear algebra, but (a) certainly parts of Sage can be improved and (b) contributions are welcome. Others can speak more knowledgeably about this particular issue.

[Thierry Dumont]

When doing numérical computations, you must use RDF floats, except for
very special purpose. RDF is your processor's float. The default Real
Field, based on MPFR (which by default as the same precision as RDF) is
slow as it is based on MPFR, a very good library, but MPFR is a software
based library.
When you make elementary linear algebra computations (factorization
eigenvalues, ...) with RDF floats, Sage uses the Lapack library based on
OpenBlas (Basic linear algebra subroutine), which is optimized for your
computer. These libraries use all your processor capacities (like the
vectorization with avx(2) instructions).
It is extremely difficult to compute faster than Lapack + OpenBlas (it
could be a bit faster with Intel own Blas, on Intels's processors...just
a bit).
To see this just compute an LU factorization:

sage: A=random_matrix(RDF,1000,1000)
sage: %time Q=A.LU()
The cost for an LU factorization is n^3 for a matrix of size n.
On my old processor (3 Ghz, avx), it takes 53.2 ms, so factorization is
computed with a speed of about 18 Gigafloats, which is impressive.

You must use RealField with a large precision (say RealField(prec=500))
if you want to solve very ill conditioned problems (example: Hilbert
matrix Aij= 1./(1.+i+j) with size n >20), but then you have to "pay"
with a slow computation.

TINA (there is no alternative :-( ), sadly.

yours,
t.d.

> <https://stackoverflow.com/questions/72206172/what-is-the-time-complexity-of-numpy-linalg-det> and made minor modifications to the script to benchmark for Sagemath.

> numpy.pngThis is the testing with matrix(RR) with the df
> algorithm.

> df.png
> I also tested with the "hessenberg" algorithm which got worse.

> df.png
> Compared to Numpy, Sagemath spent 800 seconds calculating a
> 300*300 matrix. It seems that Sagemath still runs the
> algorithm under O(n^3) (the curves of the cubic model and
> the quartic model overlap) but with a large constant. Why
> Sagemath is so slow? Is it possible to speed up?
>

> --
> 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...@googlegroups.com
> <mailto:sage-support...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sage-support/8b81830e-d76f-4554-b4cd-7f62a47ba32fn%40googlegroups.com <https://groups.google.com/d/msgid/sage-support/8b81830e-d76f-4554-b4cd-7f62a47ba32fn%40googlegroups.com?utm_medium=email&utm_source=footer>.