.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/linalg/eigen/plot_01_cauchy_spectral_theorem.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_eigen_plot_01_cauchy_spectral_theorem.py>`
        to download the full example code.

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

.. _sphx_glr_api_gallery_linalg_eigen_plot_01_cauchy_spectral_theorem.py:


Cauchy's spectral theorem: principal axes of a quadratic form
=============================================================

Cauchy's 1829 memoir showed that a real symmetric matrix has only real
eigenvalues and an orthonormal basis of eigenvectors, so
``A = V diag(lambda) V^T``. Geometrically, the eigenvectors are the
principal axes of the quadric ``x^T A x = 1``. A non-symmetric matrix
enjoys neither property.

.. GENERATED FROM PYTHON SOURCE LINES 13-18

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.linalg import eigen_general, eigen_symmetric








.. GENERATED FROM PYTHON SOURCE LINES 19-24

Principal axes of an ellipse
----------------------------
The quadratic form ``5 x^2 + 4 x y + 2 y^2 = 1`` is an ellipse whose axes
point along the eigenvectors of its symmetric matrix, with semi-axis
lengths ``1 / sqrt(lambda)``.

.. GENERATED FROM PYTHON SOURCE LINES 24-45

.. code-block:: Python


    A = np.array([[5.0, 2.0], [2.0, 2.0]])
    eig = eigen_symmetric(A)
    lams, V = eig.eigenvalues, eig.eigenvectors
    print("eigenvalues:", lams)
    print("V^T V = I:", np.allclose(V.T @ V, np.eye(2)))

    theta = np.linspace(0.0, 2.0 * np.pi, 400)
    circle = np.vstack([np.cos(theta), np.sin(theta)])
    ellipse = V @ np.diag(1.0 / np.sqrt(lams)) @ circle  # x = V Lambda^{-1/2} u, |u| = 1

    fig, ax = plt.subplots(figsize=(5, 5))
    ax.plot(ellipse[0], ellipse[1], label=r"$x^T A x = 1$")
    for k in range(2):
        axis = V[:, k] / np.sqrt(lams[k])
        ax.plot([-axis[0], axis[0]], [-axis[1], axis[1]], lw=2, label=rf"axis along $v_{k + 1}$, $\lambda_{k + 1} = {lams[k]:.0f}$")
    ax.set_aspect("equal")
    ax.set_title("Eigenvectors of a symmetric matrix = principal axes")
    ax.legend(loc="lower right", fontsize=8)
    fig.tight_layout()




.. image-sg:: /api/gallery/linalg/eigen/images/sphx_glr_plot_01_cauchy_spectral_theorem_001.png
   :alt: Eigenvectors of a symmetric matrix = principal axes
   :srcset: /api/gallery/linalg/eigen/images/sphx_glr_plot_01_cauchy_spectral_theorem_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    eigenvalues: [1. 6.]
    V^T V = I: True




.. GENERATED FROM PYTHON SOURCE LINES 46-48

A larger symmetric matrix: real spectrum, orthonormal eigenvectors
------------------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 48-58

.. code-block:: Python


    rng = np.random.default_rng(2)
    M = rng.normal(size=(6, 6))
    S = M + M.T
    eig_s = eigen_symmetric(S)
    V = eig_s.eigenvectors
    print("eigenvalues (all real):", np.round(eig_s.eigenvalues, 6))
    print(f"||V^T V - I|| = {np.linalg.norm(V.T @ V - np.eye(6)):.1e}")
    print(f"||V diag(lambda) V^T - S|| = {np.linalg.norm(V @ np.diag(eig_s.eigenvalues) @ V.T - S):.1e}")





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

 .. code-block:: none

    eigenvalues (all real): [-5.874262 -3.308306  0.673425  1.955351  3.109627  3.765125]
    ||V^T V - I|| = 2.2e-15
    ||V diag(lambda) V^T - S|| = 8.1e-15




.. GENERATED FROM PYTHON SOURCE LINES 59-63

Without symmetry the theorem fails
----------------------------------
The non-symmetric part ``M`` alone has complex eigenvalues and
eigenvectors that are not orthogonal.

.. GENERATED FROM PYTHON SOURCE LINES 63-70

.. code-block:: Python


    eig_m = eigen_general(M)
    print("eigenvalues of M:", np.round(eig_m.eigenvalues, 4))
    W = eig_m.eigenvectors
    print(f"||W^H W - I|| = {np.linalg.norm(W.conj().T @ W - np.eye(6)):.2f}")

    plt.show()




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

 .. code-block:: none

    eigenvalues of M: [-2.0452+0.j     -1.2782+0.j      0.5408+1.5451j  0.5408-1.5451j
      1.2011+0.569j   1.2011-0.569j ]
    ||W^H W - I|| = 2.33





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

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


.. _sphx_glr_download_api_gallery_linalg_eigen_plot_01_cauchy_spectral_theorem.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_cauchy_spectral_theorem.ipynb <plot_01_cauchy_spectral_theorem.ipynb>`

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

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

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

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


.. only:: html

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

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