.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/linalg/svd/plot_01_svd.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_linalg_svd_plot_01_svd.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_linalg_svd_plot_01_svd.py:


Eckart-Young: the best low-rank approximation from the SVD
==========================================================

Truncating the SVD ``A = U diag(S) V^T`` to its ``k`` largest singular
values gives the best rank-``k`` approximation of ``A``: its error is
``sqrt(sigma_{k+1}^2 + ...)`` in the Frobenius norm and ``sigma_{k+1}``
in the spectral norm, and no other rank-``k`` matrix does better. The SVD
itself comes from :func:`numpy.linalg.svd` (LAPACK's stable
bidiagonalization route, not the eigendecomposition of ``A^T A``).

.. GENERATED FROM PYTHON SOURCE LINES 14-19

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.linalg import svd_decompose








.. GENERATED FROM PYTHON SOURCE LINES 20-23

A matrix with decaying singular values
--------------------------------------
A smooth "image" sampled on a 60x80 grid, plus a little noise.

.. GENERATED FROM PYTHON SOURCE LINES 23-37

.. code-block:: Python


    rng = np.random.default_rng(3)
    y, x = np.mgrid[-1:1:60j, -1:1:80j]
    A = np.exp(-(x**2 + 2 * y**2) * 3) + 0.5 * np.cos(4 * x * y) + 0.02 * rng.normal(size=x.shape)
    res = svd_decompose(A)
    U, S, Vt = res.U, res.S, res.Vt
    print("largest singular values:", np.round(S[:6], 4))
    print("||U diag(S) Vt - A|| =", f"{np.linalg.norm(U @ np.diag(S) @ Vt - A):.1e}")


    def truncate(k):
        return U[:, :k] @ np.diag(S[:k]) @ Vt[:k]






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    largest singular values: [41.0411 15.3119  2.4311  0.3185  0.2992  0.292 ]
    ||U diag(S) Vt - A|| = 3.8e-14




.. GENERATED FROM PYTHON SOURCE LINES 38-40

Truncation error equals the discarded singular values
-----------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 40-48

.. code-block:: Python


    ks = np.arange(1, 21)
    err_fro = np.array([np.linalg.norm(A - truncate(k)) for k in ks])
    err_2 = np.array([np.linalg.norm(A - truncate(k), 2) for k in ks])
    tail = np.array([np.sqrt(np.sum(S[k:] ** 2)) for k in ks])
    print("Frobenius error == sqrt(sum of discarded sigma^2):", np.allclose(err_fro, tail))
    print("spectral error  == sigma_{k+1}:                   ", np.allclose(err_2, S[ks]))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Frobenius error == sqrt(sum of discarded sigma^2): True
    spectral error  == sigma_{k+1}:                    True




.. GENERATED FROM PYTHON SOURCE LINES 49-53

No other rank-k matrix does better
----------------------------------
Compare with rank-k approximations built by projecting ``A`` onto a
random k-dimensional column space.

.. GENERATED FROM PYTHON SOURCE LINES 53-75

.. code-block:: Python


    k = 5
    competitors = []
    for _ in range(200):
        Qk, _ = np.linalg.qr(rng.normal(size=(A.shape[0], k)))
        competitors.append(np.linalg.norm(A - Qk @ (Qk.T @ A)))
    print(f"rank-{k}: SVD error {err_fro[k - 1]:.4f}, best of 200 random projections {min(competitors):.4f}")

    fig, axes = plt.subplots(1, 4, figsize=(12, 3.2))
    axes[0].semilogy(ks, err_fro, "o-", label=r"$\|A - A_k\|_F$")
    axes[0].semilogy(ks, S[ks], "s-", label=r"$\sigma_{k+1} = \|A - A_k\|_2$")
    axes[0].set_xlabel("rank k")
    axes[0].set_title("Truncation error")
    axes[0].legend(fontsize=8)
    for ax, (label, M) in zip(axes[1:], (("A", A), ("A_1", truncate(1)), (f"A_{k}", truncate(k)))):
        ax.imshow(M, cmap="viridis", vmin=A.min(), vmax=A.max())
        ax.set_title(f"${label}$")
        ax.set_xticks([])
        ax.set_yticks([])
    fig.tight_layout()

    plt.show()



.. image-sg:: /api/gallery/linalg/svd/images/sphx_glr_plot_01_svd_001.png
   :alt: Truncation error, $A$, $A_1$, $A_5$
   :srcset: /api/gallery/linalg/svd/images/sphx_glr_plot_01_svd_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    rank-5: SVD error 1.2601, best of 200 random projections 37.2171





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.066 seconds)


.. _sphx_glr_download_api_gallery_linalg_svd_plot_01_svd.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_svd.ipynb <plot_01_svd.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_01_svd.py <plot_01_svd.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_01_svd.zip <plot_01_svd.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_