आइगेनमान और आइगेनसदिश

A <- matrix(c(4, 1,
              2, 3), nrow = 2, byrow = TRUE)
eig <- eigen(A)

lambda1 <- eig$values[1]   # 5
v1      <- eig$vectors[, 1]

# Verify: Av - lambda*v should be ~0
residual <- A %*% v1 - lambda1 * v1
max(abs(residual))  # ~ 1e-16

# Check all eigenpairs at once
# A V = V diag(lambda) where V = eigenvectors matrix
V <- eig$vectors
Lambda <- diag(eig$values)
err <- A %*% V - V %*% Lambda
max(abs(err))  # near zero

# Symmetric matrix
S <- matrix(c(4, 2, 2,
              2, 3, 1,
              2, 1, 5), nrow = 3, byrow = TRUE)

# All symmetric: A = t(A)
all.equal(S, t(S))  # TRUE

# Eigenvalues are real
eig <- eigen(S)
eig$values  # All real numbers
Im(eig$values)  # All zero

# Eigenvectors are orthogonal: V'V = I
V <- eig$vectors
round(t(V) %*% V, 10)  # Identity matrix
cat('Symmetric -> real eigenvalues, orthogonal eigenvectors')

# Rotation-like matrix -> complex eigenvalues
A <- matrix(c(0, -1,
              1,  0), nrow = 2, byrow = TRUE)
# This is a 90-degree rotation matrix
eig <- eigen(A)
eig$values
# [1] 0+1i 0-1i  (pure imaginary!)

# Real part: zero -> neither grows nor decays
Re(eig$values)  # 0 0

# Imaginary part: indicates rotation
Im(eig$values)  # 1 -1

# Another example: stable spiral
B <- matrix(c(-1, -2,
               2, -1), nrow = 2, byrow = TRUE)
eigen(B)$values  # -1 +/- 2i

A <- matrix(c(4, 1,
              2, 3), nrow = 2, byrow = TRUE)
eig <- eigen(A)
V      <- eig$vectors
Lambda <- diag(eig$values)

# Reconstruct A = V %*% Lambda %*% solve(V)
A_reconstructed <- V %*% Lambda %*% solve(V)
all.equal(A, A_reconstructed)  # TRUE

# Matrix power A^3 using eigendecomposition
A_cubed_eig <- V %*% diag(eig$values^3) %*% solve(V)
A_cubed_direct <- A %*% A %*% A
all.equal(Re(A_cubed_eig), A_cubed_direct)  # TRUE

set.seed(42)
# Correlated 2D data
x1 <- rnorm(100)
x2 <- 0.8 * x1 + 0.6 * rnorm(100)
X <- cbind(x1, x2)

# Covariance matrix
C <- cov(X)
print(round(C, 3))

# Eigendecomposition of covariance matrix
eig <- eigen(C)
cat('Eigenvalues (variance explained):\n')
print(eig$values)
cat('PC1 direction:\n')
print(eig$vectors[, 1])

# Variance explained by PC1
prop_var <- eig$values[1] / sum(eig$values)
cat('PC1 explains:', round(100 * prop_var, 1), '%')

spectral_radius <- function(A) {
  max(Mod(eigen(A)$values))
}

# Convergent matrix: spectral radius < 1
A_conv <- matrix(c(0.5, 0.2,
                   0.1, 0.3), nrow = 2, byrow = TRUE)
spectral_radius(A_conv)  # < 1 -> iterations converge

# Divergent matrix: spectral radius > 1
A_div <- matrix(c(2, 0.5,
                  0.3, 1.5), nrow = 2, byrow = TRUE)
spectral_radius(A_div)  # > 1 -> iterations diverge

# For positive definite A, spectral radius = max eigenvalue
cat('Spectral radius determines iterative stability')

A <- matrix(c(5, 2,
              1, 4), nrow = 2, byrow = TRUE)
eig_vals <- eigen(A)$values

# det(A) = product of eigenvalues
det(A)
prod(eig_vals)  # Same!

# trace = sum of eigenvalues
sum(diag(A))
sum(eig_vals)  # Same!

# For numeric precision, use Re() on complex
A2 <- matrix(c(3, -1, 2, 5), nrow = 2, byrow = TRUE)
ev <- eigen(A2)$values
all.equal(det(A2), prod(Re(ev)), tolerance = 1e-10)
all.equal(sum(diag(A2)), sum(Re(ev)), tolerance = 1e-10)

# Check positive definiteness
is_positive_definite <- function(A) {
  # Symmetric check
  if (!isTRUE(all.equal(A, t(A)))) return(FALSE)
  all(eigen(A)$values > 0)
}

# Positive definite covariance matrix
S <- matrix(c(4, 2,
              2, 3), nrow = 2, byrow = TRUE)
is_positive_definite(S)  # TRUE
eigen(S)$values  # Both positive

# Not PD (one negative eigenvalue)
Q <- matrix(c(1, 3,
              3, 2), nrow = 2, byrow = TRUE)
is_positive_definite(Q)  # FALSE
eigen(Q)$values  # One negative

# Power iteration: finds dominant eigenvalue
power_iteration <- function(A, tol = 1e-10, max_iter = 1000) {
  n <- nrow(A)
  v <- rnorm(n); v <- v / sqrt(sum(v^2))  # random unit vector
  lambda_old <- 0
  for (i in seq_len(max_iter)) {
    w <- A %*% v
    lambda <- max(abs(w))
    v <- w / lambda
    if (abs(lambda - lambda_old) < tol) break
    lambda_old <- lambda
  }
  list(value = lambda, vector = v)
}

A <- matrix(c(4, 1, 2, 3), nrow = 2, byrow = TRUE)
result <- power_iteration(A)
result$value  # Should be ~5 (dominant eigenvalue)
eigen(A)$values[1]  # Compare

A <- matrix(c(6, 2, 2, 3), nrow = 2)
eig <- eigen(A)

# Key eigen operations:
eig$values    # eigenvalues
eig$vectors   # eigenvectors (columns)

# Verify Av = lambda*v
v1 <- eig$vectors[, 1]
max(abs(A %*% v1 - eig$values[1] * v1))  # ~0

# Properties
all.equal(det(A), prod(eig$values))       # TRUE
all.equal(sum(diag(A)), sum(eig$values))  # TRUE

# Reconstruct A = V Lambda V^{-1}
V <- eig$vectors
all.equal(A, Re(V %*% diag(eig$values) %*% solve(V)))