Python version of matlab code is very slow

Question

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