Complex Eigenvalue Problems — EVLCG, EVCCG, EVLCH, EVCCH, EVLRH

This example covers five IMSL-compatible routines for complex eigenvalue problems using three physical scenarios:

• Part A — 3×3 rotation matrix: eigenvalues lie on the unit circle in the complex plane; uses eigensystems.evlcg() and eigensystems.evccg().

• Part B — 2×2 complex Hermitian matrix \(H = [[2,\,1+i],\,[1-i,\,3]]\): eigenvalues are guaranteed real; uses eigensystems.evlch() and eigensystems.evcch().

• Part C — 4×4 upper Hessenberg matrix derived from a random matrix: eigenvalues computed directly without reduction; uses eigensystems.evlrh().

Example Code

"""IMSL complex eigenvalue examples: evlcg, evccg, evlch, evcch, evlrh.

Demonstrates three scenarios:
- Part A: 3x3 rotation matrix with complex eigenvalues exp(±iθ) — evlcg, evccg
- Part B: Complex Hermitian matrix H = [[2,1+i],[1-i,3]] — evlch, evcch
- Part C: Upper Hessenberg form of a random matrix — evlrh

Outputs:
- Printed tables of eigenvalues with verification
- SVG complex-plane scatter plot saved to test_output/example_imsl_complex_eigen.svg
- Console output captured to docs/_generated_text/complex_eigen_output.txt
"""
from __future__ import annotations

from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import scipy.linalg

from eigensystems import evlcg, evccg, evlch, evcch, evlrh


def _rotation_matrix(theta: float) -> np.ndarray:
    """Build a 3x3 rotation matrix about the z-axis embedded in 3D.

    The top-left 2x2 block is the standard 2D rotation; the (3,3) element
    is 1.  Expected eigenvalues: exp(+i*theta), exp(-i*theta), 1.

    Args:
        theta (float): Rotation angle in radians.

    Returns:
        np.ndarray: 3x3 real rotation matrix.
    """
    c, s = np.cos(theta), np.sin(theta)
    return np.array([
        [ c, -s, 0.0],
        [ s,  c, 0.0],
        [0.0, 0.0, 1.0],
    ])


def run_demo_complex_eigen() -> dict:
    """Run all three complex eigenvalue demonstrations.

    Args:
        None

    Returns:
        dict: Keys ``evlcg_eigs``, ``evccg_eigs``, ``evlch_eigs``,
            ``evcch_eigs``, ``evlrh_eigs``, ``plot_path``.
    """
    # ------------------------------------------------------------------
    # Part A — rotation matrix: complex eigenvalues exp(±iθ)
    # ------------------------------------------------------------------
    theta = np.pi / 4  # 45 degrees
    R = _rotation_matrix(theta)
    # Cast to complex so evlcg / evccg accept it
    R_c = R.astype(complex)

    res_evlcg = evlcg(R_c)
    res_evccg = evccg(R_c)

    print("\nIMSL EVLCG / EVCCG Example: 3×3 Rotation Matrix")
    print(f"Rotation angle θ = π/4 = {theta:.6f} rad")
    print(f"Expected eigenvalues: exp(+iθ)={np.exp(1j*theta):.4f}, "
          f"exp(-iθ)={np.exp(-1j*theta):.4f}, 1+0j")
    print("-" * 55)
    print(f"{'EVLCG eigenvalues':}")
    for ev in sorted(res_evlcg.eigenvalues, key=lambda z: z.real):
        print(f"  {ev.real:+.6f} {ev.imag:+.6f}i")
    print(f"{'EVCCG eigenvalues (with vectors)':}")
    for i, ev in enumerate(sorted(res_evccg.eigenvalues, key=lambda z: z.real)):
        print(f"  λ{i}: {ev.real:+.6f} {ev.imag:+.6f}i")
    # Verify: |eigenvalue| should equal 1 for a rotation matrix
    mags = np.abs(res_evlcg.eigenvalues)
    print(f"Magnitude of all eigenvalues (should be 1): {mags}")

    # ------------------------------------------------------------------
    # Part B — complex Hermitian matrix H = [[2,1+i],[1-i,3]]
    # ------------------------------------------------------------------
    H = np.array([[2.0 + 0j, 1.0 + 1j],
                  [1.0 - 1j, 3.0 + 0j]])

    res_evlch = evlch(H)
    res_evcch = evcch(H)

    print("\nIMSL EVLCH / EVCCH Example: 2×2 Complex Hermitian Matrix")
    print("H = [[2, 1+i], [1-i, 3]]")
    print("Expected: eigenvalues are real (Hermitian guarantee)")
    print("-" * 55)
    print(f"EVLCH eigenvalues (values only): {res_evlch.eigenvalues}")
    print(f"EVCCH eigenvalues: {res_evcch.eigenvalues}")
    print("EVCCH eigenvectors (columns):")
    print(res_evcch.eigenvectors)
    # Verify orthonormality: V^H V = I
    V = res_evcch.eigenvectors
    orth_err = np.max(np.abs(V.conj().T @ V - np.eye(2)))
    print(f"Orthonormality error ‖V^H V - I‖_max = {orth_err:.2e}")
    # Verify all eigenvalues are real
    imag_max = np.max(np.abs(res_evlch.eigenvalues.imag))
    print(f"Max imaginary part of Hermitian eigenvalues: {imag_max:.2e} (should be ~0)")

    # ------------------------------------------------------------------
    # Part C — upper Hessenberg matrix eigenvalues via evlrh
    # ------------------------------------------------------------------
    rng = np.random.default_rng(42)
    A_rand = rng.standard_normal((4, 4))
    # Reduce to upper Hessenberg form
    H_hess, _ = scipy.linalg.hessenberg(A_rand, calc_q=True)
    # Zero out below-subdiagonal entries to guarantee Hessenberg form
    for i in range(2, 4):
        H_hess[i, :i - 1] = 0.0

    res_evlrh = evlrh(H_hess)
    eigs_direct = np.linalg.eigvals(A_rand)

    print("\nIMSL EVLRH Example: Upper Hessenberg Matrix")
    print("4×4 random matrix reduced to Hessenberg form")
    print("-" * 55)
    print(f"EVLRH eigenvalues of Hessenberg matrix:")
    for ev in sorted(res_evlrh.eigenvalues, key=lambda z: z.real):
        print(f"  {ev.real:+.6f} {ev.imag:+.6f}i")
    print(f"Direct eigvals of original matrix (for reference):")
    for ev in sorted(eigs_direct, key=lambda z: z.real):
        print(f"  {ev.real:+.6f} {ev.imag:+.6f}i")

    # ------------------------------------------------------------------
    # Plot: complex plane — all eigenvalue sets
    # ------------------------------------------------------------------
    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_complex_eigen.svg"

    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    # Subplot A — rotation matrix eigenvalues
    ax = axes[0]
    circle = plt.Circle((0, 0), 1, fill=False, linestyle="--",
                         color="gray", linewidth=0.8)
    ax.add_patch(circle)
    evs_a = res_evlcg.eigenvalues
    ax.scatter(evs_a.real, evs_a.imag, color="#be123c", s=80, zorder=5,
               label="EVLCG eigenvalues")
    ax.axhline(0, color="k", linewidth=0.4)
    ax.axvline(0, color="k", linewidth=0.4)
    ax.set_xlim(-1.5, 1.5)
    ax.set_ylim(-1.5, 1.5)
    ax.set_aspect("equal")
    ax.set_xlabel("Re(λ)")
    ax.set_ylabel("Im(λ)")
    ax.set_title("Part A — Rotation Matrix\n(EVLCG)")
    ax.legend(fontsize=8)

    # Subplot B — Hermitian eigenvalues (on real axis)
    ax = axes[1]
    evs_b = res_evcch.eigenvalues
    ax.scatter(evs_b.real, np.zeros_like(evs_b.real), color="#0369a1",
               s=120, zorder=5, label="EVCCH eigenvalues")
    ax.axhline(0, color="k", linewidth=0.4)
    ax.axvline(0, color="k", linewidth=0.4)
    for j, ev in enumerate(evs_b):
        ax.annotate(f"λ{j+1}={ev.real:.3f}", (ev.real, 0),
                    textcoords="offset points", xytext=(0, 12), fontsize=8,
                    ha="center")
    ax.set_xlim(0, 5)
    ax.set_ylim(-1, 1)
    ax.set_xlabel("Re(λ)  [Im(λ) ≈ 0]")
    ax.set_ylabel("Im(λ)")
    ax.set_title("Part B — Hermitian Matrix\n(EVLCH / EVCCH)")
    ax.legend(fontsize=8)

    # Subplot C — Hessenberg eigenvalues
    ax = axes[2]
    evs_c = res_evlrh.eigenvalues
    ax.scatter(evs_c.real, evs_c.imag, color="#059669", s=80, zorder=5,
               marker="D", label="EVLRH eigenvalues")
    ax.axhline(0, color="k", linewidth=0.4)
    ax.axvline(0, color="k", linewidth=0.4)
    ax.set_xlabel("Re(λ)")
    ax.set_ylabel("Im(λ)")
    ax.set_title("Part C — Hessenberg Matrix\n(EVLRH)")
    ax.legend(fontsize=8)

    fig.suptitle("Complex Eigenvalue Problems — IMSL Routines", fontsize=13,
                 fontweight="bold")
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)
    print(f"\nPlot saved: {plot_path}")

    return {
        "evlcg_eigs": res_evlcg.eigenvalues,
        "evccg_eigs": res_evccg.eigenvalues,
        "evlch_eigs": res_evlch.eigenvalues,
        "evcch_eigs": res_evcch.eigenvalues,
        "evlrh_eigs": res_evlrh.eigenvalues,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_complex_eigen()

Plot Output

Left: Rotation matrix eigenvalues on the unit circle (EVLCG). Centre: Hermitian eigenvalues on the real axis (EVCCH). Right: Hessenberg matrix eigenvalues in the complex plane (EVLRH).

Console Output

IMSL EVLCG / EVCCG Example: 3×3 Rotation Matrix
Rotation angle θ = π/4 = 0.785398 rad
Expected eigenvalues: exp(+i╬╕)=0.7071+0.7071j, exp(-i╬╕)=0.7071-0.7071j, 1+0j
-------------------------------------------------------
EVLCG eigenvalues
  +0.707107 +0.707107i
  +0.707107 -0.707107i
  +1.000000 +0.000000i
EVCCG eigenvalues (with vectors)
  ╬╗0: +0.707107 +0.707107i
  ╬╗1: +0.707107 -0.707107i
  ╬╗2: +1.000000 +0.000000i
Magnitude of all eigenvalues (should be 1): [1. 1. 1.]

IMSL EVLCH / EVCCH Example: 2×2 Complex Hermitian Matrix
H = [[2, 1+i], [1-i, 3]]
Expected: eigenvalues are real (Hermitian guarantee)
-------------------------------------------------------
EVLCH eigenvalues (values only): [1. 4.]
EVCCH eigenvalues: [1. 4.]
EVCCH eigenvectors (columns):
[[-0.81649658+0.j         -0.57735027+0.j        ]
 [ 0.40824829-0.40824829j -0.57735027+0.57735027j]]
Orthonormality error ΓÇûV^H V - IΓÇû_max = 4.44e-16
Max imaginary part of Hermitian eigenvalues: 0.00e+00 (should be ~0)

IMSL EVLRH Example: Upper Hessenberg Matrix
4×4 random matrix reduced to Hessenberg form
-------------------------------------------------------
EVLRH eigenvalues of Hessenberg matrix:
  -2.177794 +0.000000i
  -0.565754 +0.000000i
  +0.321431 +0.000000i
  +1.444761 +0.000000i
Direct eigvals of original matrix (for reference):
  -2.177794 +0.000000i
  -0.565754 +0.000000i
  +0.321431 +0.000000i
  +1.444761 +0.000000i

Plot saved: test_output\example_imsl_complex_eigen.svg