Hi, first time posting here. I need to compute the pseudo-inverse of large matrix (2504, 2502) as a part of an algorithm I am developing in Python (solving Ax = b). In Python, the pseudo-inverse is a clear bottleneck in terms of computation time, so I decided to look at Julia as an alternative, to check if it could improve the times.

However, when benchmarking with Python I am actually seeing slower times when using Julia. I read that LinearAlgebra uses openBLAS so I downloaded MKL.jl package, which enforces MKL usage instead of openBLAS, but even so, I seem to get slower computation times than with python/numpy.

Python:

```
import numpy as np
A = np.random.rand(2504, 2502)
```

```
%%timeit -r 5 -n 10
A_pinv = np.linalg.pinv(A)
```

leads to:

```
6.59 s Ā± 630 ms per loop (mean Ā± std. dev. of 5 runs, 10 loops each)
```

And Julia:

```
using LinearAlgebra
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 800 # so it runs all samples and evals
# BLAS.vendor() outputs :mkl
A = rand(2504, 2502)
```

```
@benchmark LinearAlgebra.pinv(x) setup=(x=$A) samples=10 evals=5
```

leads to:

```
BenchmarkTools.Trial:
memory estimate: 382.53 MiB
allocs estimate: 29
--------------
minimum time: 7.398 s (0.61% GC)
median time: 7.974 s (0.73% GC)
mean time: 8.687 s (0.69% GC)
maximum time: 10.871 s (0.65% GC)
--------------
samples: 10
evals/sample: 5
```

I also ran this more than once to make sure it wasnāt including āloading timeā for the benchmarking tools.

Any clues why this is so slow? Should I be using another package or function?