.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/linalg/stability/plot_01_condition_number.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_linalg_stability_plot_01_condition_number.py: Turing's condition number: how much a linear system amplifies error =================================================================== ``kappa_2(A) = sigma_max / sigma_min`` bounds how much a relative perturbation of the data can grow in the solution of ``A x = b``, whichever algorithm solves it: expect to lose about ``log10(kappa)`` digits. The same number explains why least squares through the normal equations ``A^T A x = A^T b`` is fragile -- ``kappa(A^T A) = kappa(A)^2`` -- while QR works with ``kappa(A)`` itself. .. GENERATED FROM PYTHON SOURCE LINES 14-22 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from scipy.linalg import hilbert from mathematicskit.linalg import condition_number_2norm, least_squares_normal_equations, least_squares_qr, lu_solve_system, svd_decompose eps = np.finfo(float).eps .. GENERATED FROM PYTHON SOURCE LINES 23-25 The condition number is a ratio of singular values -------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 25-31 .. code-block:: Python H = hilbert(6) S = svd_decompose(H).S print(f"sigma_max / sigma_min = {S.max() / S.min():.4e}") print(f"condition_number_2norm = {condition_number_2norm(H):.4e}") .. rst-class:: sphx-glr-script-out .. code-block:: none sigma_max / sigma_min = 1.4951e+07 condition_number_2norm = 1.4951e+07 .. GENERATED FROM PYTHON SOURCE LINES 32-38 Error amplification: the bound is attained ------------------------------------------ Take ``x`` along the top right singular vector and perturb ``b`` by a tiny relative amount along the worst direction (the left singular vector of ``sigma_min``): the relative change in ``x`` is exactly ``kappa`` times larger. .. GENERATED FROM PYTHON SOURCE LINES 38-47 .. code-block:: Python res = svd_decompose(H) x = res.Vt[np.argmax(res.S)] b = H @ x db = 1e-10 * np.linalg.norm(b) * res.U[:, np.argmin(res.S)] dx = lu_solve_system(H, b + db) - x amplification = (np.linalg.norm(dx) / np.linalg.norm(x)) / (np.linalg.norm(db) / np.linalg.norm(b)) print(f"relative error amplification = {amplification:.4e} (kappa = {condition_number_2norm(H):.4e})") .. rst-class:: sphx-glr-script-out .. code-block:: none relative error amplification = 1.4951e+07 (kappa = 1.4951e+07) .. GENERATED FROM PYTHON SOURCE LINES 48-52 Digits lost on Hilbert matrices ------------------------------- Solving ``H_n x = H_n 1`` in double precision: the forward error tracks ``kappa(H_n) * eps``. .. GENERATED FROM PYTHON SOURCE LINES 52-62 .. code-block:: Python ns = np.arange(2, 13) kappas, errors = [], [] for n in ns: Hn = hilbert(n) xn = lu_solve_system(Hn, Hn @ np.ones(n)) kappas.append(condition_number_2norm(Hn)) errors.append(np.linalg.norm(xn - 1.0) / np.sqrt(n)) print(f"n = {n:2d}: kappa = {kappas[-1]:.2e}, relative error = {errors[-1]:.2e}") .. rst-class:: sphx-glr-script-out .. code-block:: none n = 2: kappa = 1.93e+01, relative error = 4.23e-16 n = 3: kappa = 5.24e+02, relative error = 7.61e-16 n = 4: kappa = 1.55e+04, relative error = 1.21e-13 n = 5: kappa = 4.77e+05, relative error = 1.18e-11 n = 6: kappa = 1.50e+07, relative error = 8.63e-12 n = 7: kappa = 4.75e+08, relative error = 4.01e-09 n = 8: kappa = 1.53e+10, relative error = 4.88e-07 n = 9: kappa = 4.93e+11, relative error = 2.31e-06 n = 10: kappa = 1.60e+13, relative error = 4.55e-04 n = 11: kappa = 5.22e+14, relative error = 4.16e-03 n = 12: kappa = 1.72e+16, relative error = 2.61e-01 .. GENERATED FROM PYTHON SOURCE LINES 63-67 Least squares: kappa for QR, kappa squared for the normal equations ------------------------------------------------------------------- Build 40x5 design matrices with prescribed condition number and compare the two solvers' forward errors. .. GENERATED FROM PYTHON SOURCE LINES 67-106 .. code-block:: Python rng = np.random.default_rng(6) U, _ = np.linalg.qr(rng.normal(size=(40, 5))) V, _ = np.linalg.qr(rng.normal(size=(5, 5))) x_true = np.array([1.0, -2.0, 0.5, 3.0, -1.5]) targets = np.logspace(1, 8, 15) err_normal, err_qr = [], [] for kappa in targets: A = U @ np.diag(np.logspace(0, -np.log10(kappa), 5)) @ V.T b = A @ x_true for solver, store in ((least_squares_normal_equations, err_normal), (least_squares_qr, err_qr)): coef = solver(A, b).coefficients store.append(np.linalg.norm(coef - x_true) / np.linalg.norm(x_true)) A = U @ np.diag(np.logspace(0, -4, 5)) @ V.T print(f"kappa(A) = {condition_number_2norm(A):.3e}, kappa(A^T A) = {condition_number_2norm(A.T @ A):.3e}") fig, axes = plt.subplots(1, 2, figsize=(10, 4)) ax = axes[0] ax.loglog(kappas, errors, "o-", label="LU solve on $H_n$") ax.loglog(kappas, np.array(kappas) * eps, "k--", label=r"$\kappa\,\varepsilon$") ax.set_xlabel(r"$\kappa_2(H_n)$") ax.set_ylabel("relative forward error") ax.set_title("Hilbert systems lose log10(kappa) digits") ax.legend() ax = axes[1] ax.loglog(targets, err_normal, "o-", label="normal equations") ax.loglog(targets, err_qr, "s-", label="QR") ax.loglog(targets, targets * eps, "k--", lw=1, label=r"$\kappa\,\varepsilon$") ax.loglog(targets, targets**2 * eps, "k:", lw=1, label=r"$\kappa^2\varepsilon$") ax.set_ylim(1e-17, 10) ax.set_xlabel(r"$\kappa_2(A)$") ax.set_ylabel("relative forward error") ax.set_title("Least squares: QR vs. normal equations") ax.legend(fontsize=8) fig.tight_layout() plt.show() .. image-sg:: /api/gallery/linalg/stability/images/sphx_glr_plot_01_condition_number_001.png :alt: Hilbert systems lose log10(kappa) digits, Least squares: QR vs. normal equations :srcset: /api/gallery/linalg/stability/images/sphx_glr_plot_01_condition_number_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none kappa(A) = 1.000e+04, kappa(A^T A) = 1.000e+08 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.141 seconds) .. _sphx_glr_download_api_gallery_linalg_stability_plot_01_condition_number.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_01_condition_number.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_01_condition_number.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_01_condition_number.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_