0

I have the following Matlab code (my production code array size is of order of 8766 x 8766):

% Define z1 and z2
z1 = rand(6, 6);   % Random 6x6 matrix - real life size 8766 x 8766
z2 = rand(6, 7);   % Random 6x7 matrix - real life size 8766 x 8767

% Solve for z3
z3 = z1 \ z2; % real life size 8766 x 8767

% Display z3
disp(z3);

The equivalent of this in python is:

import numpy as np
# Define z1 and z2
z1 = np.random.rand(6, 6)   # Random 6x6 matrix
z2 = np.random.rand(6, 7)   # Random 6x7 matrix
# Solve for z3 (matrix division)
z3 = np.linalg.solve(z1.T @ z1, z1.T @ z2)

However, python version doesn't work in case it is singular matrix. The Matlab version works even if matrix is singular.

And the following alternate solution is tediously slow:

z3 = np.linalg.pinv(z1) @ z2

So, what is the solution for me?

Amazing that in that in 2023, Python still hasn't beat Matlab on that.

Edit: scipy.linalg.lstsq is slower than methods already used

Pranav Hosangadi
  • 23,755
  • 7
  • 44
  • 70
Zanam
  • 4,607
  • 13
  • 67
  • 143
  • Did you try [scipy.linalg.lstsq](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html#scipy.linalg.lstsq)? – Cris Luengo Jun 29 '23 at 19:10
  • 2
    `z3 = z1 / z2` is wrong, it would be the solution to `z3 * z2 = z1`. Maybe you mean to do `z2 \ z1`? Can you show your timing code? With these sizes, the solution is found instantaneously. Maybe you're using larger arrays? – Cris Luengo Jun 30 '23 at 07:34
  • @Cris yes my real life array is much large and I fixed the Z2/Z1 situation. – Zanam Jun 30 '23 at 14:37
  • 1
    Your Python code solves for `z3` in `z2 * z3 = z1`. Your updated MATLAB code solves for `z3` in `z3 * z1 = z2`. If you can’t match the operations, why would we trust your timings? (I don’t doubt MATLAB is faster, it’s really good at picking the right algorithms) Please post the code you used for timing, it’ll clarify what you are trying to achieve and help us test possible solutions. – Cris Luengo Jun 30 '23 at 14:56
  • 1
    Can you time them and add the actual, measured speed differences? That would be interesting data. – Gabriel Staples Jun 30 '23 at 15:03
  • See [here](https://stackoverflow.com/a/55483480/7328782) for some pointers. With your larger array sizes, I see 21s in MATLAB, and 69s in Python for `z3, _, _, _ = scipy.linalg.lstsq(z1, z2, lapack_driver='gelsy')`. Yes, it's slower, but it's not awfully slow. Using a faster LAPACK implementation in Python will likely make it as fast as in MATLAB. Google should show you how to switch out the LAPACK library. – Cris Luengo Jun 30 '23 at 15:46
  • @Cris can you please elaborate a bit more on LAPACK solution – Zanam Jun 30 '23 at 15:52
  • https://stackoverflow.com/questions/9000164/how-to-check-blas-lapack-linkage-in-numpy-and-scipy – Pranav Hosangadi Jun 30 '23 at 15:54
  • @CrisLuengo the solution feels awfully slow since in my application I have to run it for 1000s of different instances. I am still surprised that it is not a default implementation in Python. I feel like I am going to reinvent the wheel. – Zanam Jul 01 '23 at 12:39

0 Answers0