# 单位圆经R变换成椭圆
theta <- seq(0, 2*pi, length=100)
circle <- cbind(cos(theta), sin(theta))

# 假设 R = [[4, 2], [2, 3]]
R <- matrix(c(4, 2, 2, 3), nrow=2)
ellipse <- circle %*% R

# 绘图展示
plot(circle, type="l", asp=1, main="单位圆→椭圆变换")
lines(ellipse, col="red", lwd=2)

# 生成对称矩阵
set.seed(123)
m <- 5
R <- matrix(runif(m^2, 0.1, 0.9), m, m)
R <- (R + t(R)) / 2  # 强制对称
diag(R) <- 1         # 对角线设为1（LD矩阵）

# 特征值分解
eig <- eigen(R, symmetric = TRUE)

# 提取结果
U <- eig$vectors      # 特征向量矩阵
values <- eig$values  # 特征值向量
Lambda <- diag(values) # 特征值对角矩阵

# 验证分解正确性
R_recon <- U %*% Lambda %*% t(U)
all.equal(R, R_recon, tolerance = 1e-10)  # 应返回TRUE

# 只保留前k个主成分
k <- 2  # 保留的主成分数量

# 选择前k个特征值和特征向量
U_k <- U[, 1:k]                # m × k
values_k <- values[1:k]        # 长度k
Lambda_k <- diag(values_k)     # k × k

# 低秩近似重建
R_approx <- U_k %*% Lambda_k %*% t(U_k)  # m × m，但秩为k

# LD矩阵的特征值分解
analyze_ld_structure <- function(R, eigen_ratio = 0.95) {
    eig <- eigen(R, symmetric = TRUE)
    values <- eig$values
    U <- eig$vectors
    
    # 计算累积方差解释比例
    total_var <- sum(values)
    cum_prop <- cumsum(values) / total_var
    
    # 确定保留的主成分数
    k <- which(cum_prop >= eigen_ratio)[1]
    if (is.na(k)) k <- length(values)
    
    # 输出结果
    cat("总方差:", total_var, "\n")
    cat("特征值分布:\n")
    print(head(values, 10))
    cat("\n保留主成分数:", k, "\n")
    cat("解释方差比例:", cum_prop[k], "\n")
    
    return(list(U = U, values = values, k = k))
}

# 技巧1：计算 R × beta（避免显式构造R）
compute_R_times_beta <- function(beta, U, values) {
    # beta: m×1向量
    # U: m×k特征向量矩阵
    # values: k个特征值
    
    # 高效计算：R × beta = U × diag(values) × U^T × beta
    tmp <- crossprod(U, beta)  # = U^T × beta，k×1矩阵
    tmp <- tmp[, 1] * values   # 转换为向量并乘以特征值
    result <- U %*% tmp        # = U × (diag(values) × U^T × beta)
    
    return(result)
}

# 技巧2：计算 R^{-1} × beta（避免求逆）
compute_R_inv_times_beta <- function(beta, U, values, epsilon = 1e-8) {
    # 过滤小特征值，避免数值问题
    keep <- values > epsilon
    U_filtered <- U[, keep]
    values_filtered <- values[keep]
    
    # 计算：R^{-1} × beta = U × diag(1/values) × U^T × beta
    tmp <- crossprod(U_filtered, beta)
    tmp <- tmp[, 1] / values_filtered  # 关键：除以特征值！
    result <- U_filtered %*% tmp
    
    return(result)
}

# 确保矩阵正定
ensure_positive_definite <- function(R, eig_tol = 1e-8) {
    eig <- eigen(R, symmetric = TRUE)
    values <- eig$values
    U <- eig$vectors
    
    # 检查最小特征值
    min_eig <- min(values)
    if (min_eig < eig_tol) {
        cat(sprintf("警告：最小特征值 = %.2e，进行校正\n", min_eig))
        
        # 方法1：添加岭惩罚
        lambda <- abs(min_eig) + 1e-6
        R_corrected <- R + diag(lambda, nrow(R))
        
        # 重新标准化（如果是相关矩阵）
        s <- 1 / sqrt(diag(R_corrected))
        R_corrected <- diag(s) %*% R_corrected %*% diag(s)
        
        return(R_corrected)
    }
    
    return(R)
}

# 使用Matrix包的nearPD函数
if (!require(Matrix)) install.packages("Matrix")
library(Matrix)

correct_with_nearpd <- function(R) {
    # 寻找最近的正定矩阵
    R_pd <- nearPD(R, 
                   corr = TRUE,      # 保持相关矩阵结构
                   keepDiag = TRUE,  # 保持对角线
                   do2eigen = TRUE,  # 确保正定
                   eig.tol = 1e-6,   # 特征值容忍度
                   conv.tol = 1e-7,  # 收敛容忍度
                   maxit = 100)      # 最大迭代次数
    
    return(as.matrix(R_pd$mat))
}