A \\in \\mathbb{R}^{n,p}<\/span><\/code> and a vector b \\in \\mathbb{R}^{n}<\/span><\/code> and one is interested in solution vectors x \\in \\mathbb{R}^{p}<\/span><\/code> fulfilling<\/p>\n\n\n\nx^\\star = \\argmin_x ||Ax-b||_2^2 \\,.<\/pre><\/div>\n\n\n\nYou can read more about least squares in our earlier post Least Squares Minimal Norm Solution<\/a>. There are many possible ways to tackle this problem and many algorithms are available. One standard solver is LSQR<\/a> with the following properties:<\/p>\n\n\n\n- Iterative solver, which terminates when some stopping criteria are smaller than a user specified tolerance.<\/li>
- It only uses matrix-vector products. Therefore, it is suitable for large and sparse matrices A<\/span>.<\/li>
- It effectively solves the normal equations A^T A = A^Tb<\/span> based on a bidiagonalization procedure of Golub and Kahan (so never really calculating A^T A<\/span>).<\/li>
- It is, unfortunately, susceptible to ill-conditioned matrices A<\/span>, which we demonstrated in our earlier post.<\/li><\/ul>\n\n\n\n
Wait, what is an ill-conditioned matrix? This is most easily explained with the help of the singular value decomposition<\/a> (SVD). Any real valued matrix permits a decomposition into three parts:<\/p>\n\n\n\nA = U S V^T \\,.<\/pre><\/div>\n\n\n\nU<\/span> and V<\/span> are orthogonal matrices, but not of further interest to us. The matrix S<\/span> on the other side is (rectangular) diagonal with only non-negative entries s_i = S_{ii} \\geq 0<\/span>. A matrix A<\/span> is said to be ill-conditioned if it has a large condition number, which can be defined as the ratio of largest and smallest singular value, \\mathrm{cond}(A) = \\frac{\\max_i s_i}{\\min_i s_i} = \\frac{s_{\\mathrm{max}}}{s_{\\mathrm{min}}}<\/span>. Very often, large condition numbers make numerical computations difficult or less precise.<\/p>\n\n\n\nBenchmarking LSQR<\/h2>\n\n\n\n
One day, I decided to benchmark the computation time of least squares solvers provided by scipy<\/a>, in particular LSQR. I wanted results for different sizes n, p<\/span> of the matrix dimensions. So I needed to somehow generate different matrices A<\/span>. There are a myriad ways to do that. Naive as I was, I did the most simple thing and used standard Normal (Gaussian) distributed random matrices A_{ij} \\sim \\mathcal{N}(0, 1)<\/span> and ran benchmarks on those. Let’s see how that looks in Python.<\/p>\n\n\n\nfrom collections import OrderedDict\nfrom functools import partial\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.linalg import lstsq\nfrom scipy.sparse.linalg import lsqr\nimport seaborn as sns\nfrom neurtu import Benchmark, delayed\n\nplt.ion()\n\np_list = [100, 500]\nrng = np.random.default_rng(42)\nX = rng.standard_normal(max(p_list) ** 2 * 2)\ny = rng.standard_normal(max(p_list) * 2)\n\ndef bench_cases(): \n for p in p_list:\n for n in np.arange(0.1, 2.1, 0.1):\n n = int(p*n)\n A = X[:n*p].reshape(n, p)\n b = y[:n]\n for solver in ['lstsq', 'lsqr']:\n tags = OrderedDict(n=n, p=p, solver=solver)\n if solver == 'lstsq':\n solve = delayed(lstsq, tags=tags)\n elif solver == 'lsqr':\n solve = delayed(\n partial(\n lsqr, atol=1e-10, btol=1e-10, iter_lim=1e6\n ),\n tags=tags)\n\n yield solve(A, b)\n\nbench = Benchmark(wall_time=True)\ndf = bench(bench_cases())\n\ng = sns.relplot(x='n', y='wall_time', hue='solver', col='p',\n kind='line', facet_kws={'sharex': False, 'sharey': False},\n data=df.reset_index(), marker='o')\ng.set_titles("p = {col_name}")\ng.set_axis_labels("n", "Wall time (s)")\ng.set(xscale="linear", yscale="log")\n\nplt.subplots_adjust(top=0.9)\ng.fig.suptitle('Benchmark LSQR vs LSTSQ')\nfor ax in g.axes.flatten():\n ax.tick_params(labelbottom=True)<\/pre><\/div>\n\n\n\n![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2022\/05\/fig_benchmark-1024x475.png\")
Timing benchmark of LSQR and standard LAPACK LSTSQ solvers for different sizes of n on the x-axis and p=100 left, p=500 right.<\/figcaption><\/figure>\n\n\n\nThe left plot already looks a bit suspicious around n=p<\/span>. But what is happening on the right side? Where does this spike of LSQR come from? And why does the standard least squares solver, SVD-based lstsq, not show this spike?<\/p>\n\n\n\nWhen I saw these results, I thought something might be wrong with LSQR and opened an issue on the scipy github repository, see <\/a>https:\/\/github.com\/scipy\/scipy\/issues\/11204<\/a>. The community there is really fantastic. Brett Naul pointed me to ….<\/p>\n\n\n\nThe Marchenko\u2013Pastur Distribution<\/h2>\n\n\n\n
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2022\/05\/empirical_eigenvalue_distribution-1024x768.png\")
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2022\/05\/marchenko_pastur_distribution-1024x768.png\")
![\"\"](\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2022\/05\/condition_number-1024x768.png\")