<p><strong>Give a try to <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html#scipy.linalg.lstsq" rel="nofollow noreferrer"><code>scipy.linalg.lstsq()</code></a> using <code>lapack_driver='gelsy'</code>!</strong></p>

<p>Let's review the different routines for solving linear least square and the approches:</p>

<ul>
<li><p><a href="https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html" rel="nofollow noreferrer"><code>numpy.linalg.lstsq()</code></a> wraps LAPACK's <a href="http://www.netlib.org/lapack/explore-html/db/d6a/dgelsd_8f_source.html" rel="nofollow noreferrer"><code>xGELSD()</code></a>, as shown in <a href="https://github.com/numpy/numpy/blob/v1.16.1/numpy/linalg/umath_linalg.c.src" rel="nofollow noreferrer">umath_linalg.c.src</a> on line 2841+. This routine reduces the matrix V to bidiagonal form using a devide and conquer strategy and compute the SVD of that bidiagonal matrix.</p></li>
<li><p>scipy's <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html#scipy.linalg.lstsq" rel="nofollow noreferrer"><code>scipy.linalg.lstsq()</code></a> wraps LAPACK's <code>xGELSD()</code>, <a href="http://www.netlib.org/lapack/explore-html/d6/d4b/dgelsy_8f_source.html" rel="nofollow noreferrer"><code>xGELSY()</code></a> and <a href="http://www.netlib.org/lapack/explore-html/d9/d4e/dgelss_8f_source.html" rel="nofollow noreferrer"><code>xGELSS()</code></a>: the argument <code>lapack_driver</code> can be modifed to switch from one to another. According to the benchmark of LAPACK, <code>xGELSS()</code> is slower than <code>xGELSD()</code> and <code>xGELSY()</code> is about 5 faster than <code>xGELSD()</code>. <code>xGELSY()</code> makes use of a QR factorization of V with column pivoting. And the good news is that this switch was <a href="https://docs.scipy.org/doc/scipy-1.1.0/reference/generated/scipy.linalg.lstsq.html#scipy.linalg.lstsq" rel="nofollow noreferrer">already available in scipy 1.1.0</a>! </p></li>
<li>LAPACK's <a href="http://www.netlib.org/lapack/explore-html/d8/dde/dgels_8f_source.html" rel="nofollow noreferrer"><code>xGELS()</code></a> makes use of the QR decomposition of the matrix V, but it assumes that this matrix has full rank. According to the benchmark of LAPACK, on can expect <code>dgels()</code> to be about 5 time faster than <code>dgelsd()</code>, but it also more vulnerable to the condition number of the matrix and may become innacurate. See details and further references in <a href="https://stackoverflow.com/questions/41637108/the-difference-between-c-lapack-sgels-and-python-numpy-lstsq-results">The difference between C++ (LAPACK, sgels) and Python (Numpy, lstsq) results</a> . xGELS() is available in the <a href="https://docs.scipy.org/doc/scipy/reference/linalg.cython_lapack.html#module-scipy.linalg.cython_lapack" rel="nofollow noreferrer">cython-lapack interface of scipy</a>.  </li>
</ul>

<p>While very tempting, computing and using <code>V^T·V</code> to solve the normal equation is likely not the way to go. Indeed, the precision is endangered by the condition number of that matrix, about <a href="https://www.cs.ubc.ca/~rbridson/courses/542g-fall-2008/notes-oct1.pdf" rel="nofollow noreferrer">the square of the condition number of the matrix V</a>. Since <a href="https://arxiv.org/abs/1504.02118" rel="nofollow noreferrer">Vandermonde matrices tend to be badly ill-conditioned, except for the matrices of the discrete Fourier transform</a>, it can become hazardous... <strong>Finally, you may even keep using <code>xGELSD()</code> to avoid problems related to conditionning. If you switch to <code>xGELSY()</code>, consider <a href="http://www.netlib.org/lapack/lug/node82.html" rel="nofollow noreferrer">estimating the error</a>.</strong></p>
Fastest way of solving linear least squares
<p>In <a href="https://math.stackexchange.com/a/2233298/340174">https://math.stackexchange.com/a/2233298/340174</a> it is mentioned that solving a linear equation "M·x = b" (matrix M is square) is slow if done via LU decomposition (but even slower using QR decomposition). Now I noticed that <code>numpy.linalg.solve</code> is in fact using LU decomposition. In truth, I want to solve "V·x = b" for a non-squared Vandermonde design matrix V for the least squares. I don't want regularization. I see multiple approaches:</p>

<ol>
<li>Solve "V·x = b" with <code>numpy.linalg.lstsq</code>, which uses Fortran "xGELSD" based on SVD. The SVD should be even slower than LU decomposition, but I don't need to calculate "(V^T·V)".</li>
<li>Solve "(V^T·V)·x = (V^T·b)" with <code>numpy.linalg.solve</code>, which uses LU decomposition.</li>
<li>Solve "A·x = b" with <code>numpy.linalg.solve</code>, which uses LU decomposition, but calculating "A=xV^T·V" directly according to <a href="https://math.stackexchange.com/a/3155891/340174">https://math.stackexchange.com/a/3155891/340174</a></li>
</ol>

<p>Alaternatively I could use the newest <code>solve</code> from scipy (<a href="https://docs.scipy.org/doc/scipy-1.2.1/reference/generated/scipy.linalg.solve.html" rel="nofollow noreferrer">https://docs.scipy.org/doc/scipy-1.2.1/reference/generated/scipy.linalg.solve.html</a>) which can use diagonal pivoting for the symmetric matrix "A" (which is faster than using LU decomposition, I guess), but my scipy is stuck on 1.1.0, so I don't have access to that.</p>

<p>From <a href="https://stackoverflow.com/a/45535523/4533188">https://stackoverflow.com/a/45535523/4533188</a> it seems that <code>solve</code> is faster than <code>lstsq</code>, including calculating "V^T·V", but when I tried it, <code>lstsq</code> was faster. Maybe I am doing something wrong?</p>

<p><strong>What is the fastest way of solving my linear problem?</strong></p>

<hr>

<h3>No real options</h3>

<ul>
<li><code>statsmodels.regression.linear_model.OLS.fit</code> is uses either Moore-Penrose pseudoinverse or QR-factorization + <code>np.linalg.inv</code> + <code>np.linalg.svd</code> + <code>numpy.linalg.solve</code>, which does not seem too efficient to me.</li>
<li><code>sklearn.linear_model.LinearRegression</code> uses scipy.linalg.lstsq.</li>
<li><code>scipy.linalg.lstsq</code> uses also xGELSD.</li>
<li>I expect calculating the inverse of "(V^T·V)" to be pretty expensive, so I discarded the direct computation of "x = (V^T·V)^-1·(V^T·b)"</li>
</ul>