I have been looking at the function numpy.linalg.lstsq and it contains all types of solver that LAPACK provides except for GELS routine that is based on QR decomposition (householder reflectors in the matter of fact as computing the QR with Q is a waste of computation time). The question is there is an implementation for this GELS routine?
Moreover, in all implementations that I saw in the web for Least squares based on QR, they compute the Q in the QR decomposition, and then use it, but it needs to inefficient as computing Q is a waste, and we need only the reflectores.
the website of LAPCAK that explain what are the type of Least squares solver: https://www.netlib.org/lapack/lug/node27.html as you can see numpy.linalg.lstsq support only 3 out of the 4.
If you can't find a high-level counterpart of a LAPACK routine you want, there's always scipy.linalg.lapack, which contains low-level wrappers for LAPACK routines. It looks like all the GELS variants are there, except the GELSX compatibility variants - there's dgels, sgels, etc.