Of course! This is a very common and important topic in scientific computing with Python. Let's break down how to compute eigenvalues and eigenvectors using the two most popular libraries: NumPy and SciPy.

The Core Concept: What are Eigenvalues and Eigenvectors?

Before we dive into the code, let's quickly recall the math. For a given square matrix A, we are looking for a special scalar value (λ, the eigenvalue) and a non-zero vector (v, the eigenvector) that satisfy the equation:

A v = λ v

This equation means that when you multiply the matrix A by the vector v, the result is simply the original vector v scaled by the factor . The vector v doesn't change its "direction" in the vector space, only its magnitude.

Using NumPy

NumPy is the fundamental package for numerical computation in Python. Its linear algebra module, numpy.linalg, is the standard go-to for many matrix operations, including eigenvalue decomposition.

The Main Function: numpy.linalg.eig

This function takes a square matrix as input and returns two arrays:

1. w: A 1D array containing the eigenvalues.
2. v: A 2D array where each column is the eigenvector corresponding to the eigenvalue in w.

Important Note on Eigenvectors: NumPy returns eigenvectors that are normalized (have a length or "norm" of 1).

Example:

Let's find the eigenvalues and eigenvectors for a simple 2x2 matrix: A = [[4, 1], [2, 3]]

import numpy as np
# Define the matrix
A = np.array([[4, 1],
              [2, 3]])
# Calculate eigenvalues and eigenvectors
# The function returns a tuple: (eigenvalues, eigenvectors)
eigenvalues, eigenvectors = np.linalg.eig(A)
print("Matrix A:")
print(A)
print("\nEigenvalues (w):")
print(eigenvalues)
print("\nEigenvectors (v) - each column is an eigenvector:")
print(eigenvectors)
# --- Verification ---
# Let's check if A @ v_i = lambda_i * v_i for the first pair
# Note: eigenvectors are columns, so we use [:, 0]
lambda_1 = eigenvalues[0]
v_1 = eigenvectors[:, 0]
# Perform the matrix-vector multiplication
result_Av = np.dot(A, v_1)
result_lambda_v = lambda_1 * v_1
print("\n--- Verification for the first pair ---")
print(f"A @ v_1 = \n{result_Av}")
print(f"lambda_1 * v_1 = \n{result_lambda_v}")
print(f"Are they equal? {np.allclose(result_Av, result_lambda_v)}")

Output:

Matrix A:
[[4 1]
 [2 3]]
Eigenvalues (w):
[5. 2.]
Eigenvectors (v) - each column is an eigenvector:
[[ 0.70710678 -0.4472136 ]
 [ 0.70710678  0.89442719]]
--- Verification for the first pair ---
A @ v_1 = 
[3.53553391 3.53553391]
lambda_1 * v_1 = 
[3.53553391 3.53553391]
Are they equal? True

As you can see, the verification confirms that our calculation is correct.

Using SciPy

SciPy is built on top of NumPy and provides more advanced functionality for scientific and technical computing. Its linear algebra module, scipy.linalg, is often preferred for more complex or specialized tasks.

The Main Function: scipy.linalg.eig

For basic eigenvalue problems, scipy.linalg.eig works almost identically to numpy.linalg.eig. The key difference is that SciPy's functions are often highly optimized and can be more robust for certain edge cases.

Example:

The code is very similar to the NumPy version.

import numpy as np
from scipy import linalg
# Define the same matrix
A = np.array([[4, 1],
              [2, 3]])
# Calculate eigenvalues and eigenvectors using SciPy
eigenvalues_scipy, eigenvectors_scipy = linalg.eig(A)
print("Matrix A (same as before):")
print(A)
print("\nEigenvalues from SciPy:")
print(eigenvalues_scipy)
print("\nEigenvectors from SciPy:")
print(eigenvectors_scipy)
# The results should be identical to NumPy's
print("\nAre NumPy and SciPy eigenvalues equal?", np.allclose(eigenvalues, eigenvalues_scipy))
print("Are NumPy and SciPy eigenvectors equal?", np.allclose(eigenvectors, eigenvectors_scipy))

Output:

Matrix A (same as before):
[[4 1]
 [2 3]]
Eigenvalues from SciPy:
[5.+0.j 2.+0.j]
Eigenvectors from SciPy:
[[ 0.70710678+0.j -0.4472136 +0.j]
 [ 0.70710678+0.j  0.89442719+0.j]]
Are NumPy and SciPy eigenvalues equal? True
Are NumPy and SciPy eigenvectors equal? True

Note that SciPy returns eigenvalues and eigenvectors with a complex data type (dtype=complex128), even for real matrices. This is because it's designed to handle complex matrices from the start, which is useful for more advanced problems.

When to Use NumPy vs. SciPy?

Recommendation:

• For most standard eigenvalue problems, numpy.linalg.eig is perfectly fine and is the standard choice.
• If you are working with symmetric or Hermitian matrices, you should use scipy.linalg.eigh. It's faster and more numerically stable for this specific case.
• If you are working with very large, sparse matrices, you should use scipy.sparse.linalg.eigsh, which can find only a few eigenvalues/vectors without computing all of them, saving a massive amount of time and memory.

Important Special Case: Hermitian/Symmetric Matrices

For a real, symmetric matrix (or a complex Hermitian matrix), all eigenvalues are real numbers. The functions numpy.linalg.eigh and scipy.linalg.eigh are designed for this. They are more efficient and stable.

# A symmetric matrix
S = np.array([[1, 2, 3],
              [2, -4, 5],
              [3, 5, 6]])
# Using SciPy's eigh (or NumPy's)
# Note: eigh takes advantage of the matrix's symmetry
eigenvalues_sym, eigenvectors_sym = linalg.eigh(S)
print("Symmetric Matrix S:")
print(S)
print("\nEigenvalues (guaranteed to be real):")
print(eigenvalues_sym)
print("\nEigenvectors:")
print(eigenvectors_sym)

Summary