基于Power Method求解矩阵全部特征值与特征向量的R语言咨询

Hey Joel, great job getting the Power Method working for the dominant eigenvalue and eigenvector—half the battle is already won! Recursion paired with matrix deflation is exactly the elegant approach you’re looking for to generalize this to all eigenvalues and eigenvectors. Here’s how to make it happen:

Core Idea #

The plan is straightforward:

1. Use your existing Power Method to find the dominant (largest magnitude) eigenvalue and eigenvector of the current matrix.
2. "Deflate" the matrix to remove the component tied to this eigenvalue, leaving a smaller matrix that holds the remaining eigenvalues.
3. Recursively repeat steps 1-2 on the deflated matrix until we’ve extracted every eigenvalue and corresponding eigenvector.

Step 1: Refine Your Power Method Function #

First, let’s turn your working Power Method into a reusable function that returns both the eigenvalue and normalized eigenvector:

power_method <- function(A, tol = 1e-6, max_iter = 1000) {
  n <- nrow(A)
  x <- runif(n)  # Start with a random initial vector
  x <- x / norm(x, "2")  # Normalize to prevent overflow
  
  for (i in 1:max_iter) {
    x_new <- A %*% x
    lambda <- t(x) %*% x_new  # Rayleigh quotient for eigenvalue estimate
    x_new <- x_new / norm(x_new, "2")
    
    # Check if we've converged
    if (norm(x_new - x, "2") < tol) {
      break
    }
    x <- x_new
  }
  
  return(list(eigenvalue = as.numeric(lambda), eigenvector = x))
}

Step 2: Recursive Deflation Function #

Next, we’ll build a recursive function that uses the Power Method and deflates the matrix each time. We’ll use Hotelling’s deflation—a stable technique that preserves symmetry (if your original matrix is symmetric) and keeps numerical errors in check:

recursive_eigen <- function(A, tol = 1e-6, max_iter = 1000) {
  n <- nrow(A)
  
  # Base case: 1x1 matrix, just return the single value and a trivial vector
  if (n == 1) {
    return(list(
      eigenvalues = as.numeric(A),
      eigenvectors = matrix(1, nrow = 1)
    ))
  }
  
  # Grab the dominant eigenvalue and eigenvector
  dominant <- power_method(A, tol, max_iter)
  lambda <- dominant$eigenvalue
  v <- dominant$eigenvector
  
  # Deflate the matrix to remove the dominant component
  deflated_A <- A - lambda * (v %*% t(v))
  
  # Recurse on the smaller deflated matrix
  recursive_result <- recursive_eigen(deflated_A, tol, max_iter)
  
  # Combine results: prepend the dominant pair to the recursive output
  combined_eigenvalues <- c(lambda, recursive_result$eigenvalues)
  combined_eigenvectors <- cbind(v, recursive_result$eigenvectors)
  
  return(list(
    eigenvalues = combined_eigenvalues,
    eigenvectors = combined_eigenvectors
  ))
}

Step 3: Test It Out #

Let’s validate this with a sample matrix where we know the expected eigenvalues:

# 3x3 matrix with known eigenvalues: 3, 2, 1
A <- matrix(c(4, -1, -1, -1, 3, -1, -1, -1, 2), nrow = 3)

# Run our recursive function
result <- recursive_eigen(A)

# Print the results
cat("Extracted Eigenvalues:\n")
print(result$eigenvalues)

cat("\nExtracted Eigenvectors (each column is a vector):\n")
print(result$eigenvectors)

# Compare with R's built-in eigen() function for sanity
cat("\nBuilt-in eigen() results for comparison:\n")
print(eigen(A))

Key Notes #

• Deflation Stability: Hotelling’s method is preferred over simpler deflation techniques because it avoids amplifying numerical errors, especially for symmetric matrices.
• Repeated Eigenvalues: If your matrix has repeated eigenvalues, you might need to add checks or adjust the deflation step—this implementation works best when eigenvalues are distinct in magnitude.
• Normalization: We normalize eigenvectors at every step to keep values manageable and ensure consistent convergence.

This approach builds directly on your existing Power Method code and uses recursion to elegantly handle matrices of any size (as long as the Power Method converges for each deflated matrix).

内容的提问来源于stack exchange，提问作者Joel H