手搓矩阵相关算法及函数

矩阵消元法

行式最简形（Row-Reduced Echelon Form，RREF）

矩阵消元法，这通常涉及到将矩阵转换为行最简形。然而，R的基础包并没有直接提供计算RREF的函数。不过，我们可以手动通过一系列的元素操作来模拟这个过程，或者使用额外的包（如pracma）来实现这一功能。

矩阵消元的简单演示：

# 定义一个3x3的矩阵
A <- matrix(c(2, 1, -1,
              -3, -1, 2,
              -2, 1, 2), byrow = TRUE, nrow = 3)

# 定义对应的常数项向量
b <- c(8, -11, -3)

# 将常数项向量转化为矩阵的一列，使其成为增广矩阵
Ab <- cbind(A, b)

# 手动执行高斯消元法
# 第一步：使用第一行消去以下各行的第一个元素
Ab[2,] <- Ab[2,] + 3/2 * Ab[1,]
Ab[3,] <- Ab[3,] + Ab[1,]

# 第二步：使用第二行消去以下各行的第二个元素
Ab[3,] <- Ab[3,] - Ab[2,] * Ab[3,2] / Ab[2,2]

# 输出部分消元后的增广矩阵
Ab

                b
[1,] 2 1.0 -1.0 8
[2,] 0 0.5  0.5 1
[3,] 0 0.0 -1.0 1

# 你可以继续以上步骤，直到矩阵变为上三角形矩阵
# 然后通过回代（back substitution）求解每个变量

在实际中，需要继续进行行操作直到矩阵变为上三角矩阵，然后通过回代求解未知数。

如果希望使用一个现成的函数来求解RREF，可以使用pracma包：

# 使用pracma包
library(pracma)

# 求解行最简形
rref(Ab)

            b
[1,] 1 0 0  2
[2,] 0 1 0  3
[3,] 0 0 1 -1

回代（back substitution）

假设我们已经有了一个经过高斯消元的上三角增广矩阵，下面是如何进行回代操作的示例：

# 假设我们有如下的上三角增广矩阵
Ab <- matrix(c(2, 1, -1, 8,
               0, -2.5, 3.5, -1.5,
               0, 0, -2, 2), byrow = TRUE, nrow = 3)

# 回代过程
n <- nrow(Ab)
x <- numeric(n)  # 创建一个长度为n的向量来存储解

# 从最后一行开始进行回代
x[n] <- Ab[n, n+1] / Ab[n, n]  # 最后一个方程的解
for (i in (n-1):1) {
  # 对于每一行，从增广部分减去已经求得的解与对应系数的乘积
  x[i] <- (Ab[i, n+1] - sum(Ab[i, (i+1):n] * x[(i+1):n])) / Ab[i, i]
}

# 打印解
x

[1]  3.9 -0.8 -1.0

在这个代码段中，首先创建了一个长度为n的向量x来存储解，其中n是矩阵Ab的行数。然后我们从最后一行开始，对每一行进行以下操作：

1. 从增广矩阵的最后一列中取出对应的常数项。

2. 从当前行中，减去已经解出的未知数与其对应系数的乘积。

3. 用上述结果除以当前行的主对角线元素（即未知数的系数），得到当前未知数的解。

在上述代码中，Ab[i, n+1] 是增广矩阵最后一列的第 i 行元素，代表常数项；Ab[i, i] 是第 i 行第 i 列的元素，即第 i 个未知数的系数；sum(Ab[i, (i+1):n] * x[(i+1):n]) 计算已知解对当前行的影响。

请注意，这个过程假设了矩阵已经是上三角形，并且没有零在主对角线上。在实际应用中，还需要检查对角线元素是否为零（即方程是否有唯一解），并且可能需要处理数值不稳定的情况。在 R 中，通常不需要手动进行这些步骤，因为可以直接使用 solve() 函数来求解整个方程组。

# 定义一个3x3的矩阵
A <- matrix(c(2, 1, -1,
              -3, -1, 2,
              -2, 1, 2), byrow = TRUE, nrow = 3)

# 定义对应的常数项向量
b <- c(8, -11, -3)

solve(A, b)

[1]  2  3 -1

LU 分解

LU分解是一种矩阵分解技术，它将一个给定的方阵分解为两个特殊的矩阵的乘积：一个下三角矩阵 \(L\) 和一个上三角矩阵 \(U\)。

这种分解可以表示为： \(A = LU\)

其中：

• \(A\) 是要分解的原始矩阵。
• \(L\) 是一个下三角矩阵，其对角线上的元素通常设为1。
• \(U\) 是一个上三角矩阵。

LU分解的目的是为了简化矩阵方程的求解过程，比如线性方程组 \(Ax = b\)。通过分解 \(A\) 为 \(LU\)，原方程可转换为 \(LUx = b\)。这可以通过以下步骤求解：

1. 首先求解 \(Ly = bZ\) 通过前代（forward substitution）。
2. 然后求解 \(Ux = y\) 通过回代（backward substitution）。

LU分解特别适用于需要多次求解不同右侧向量 \(b\) 的相同系数矩阵 \(A\) 的方程组，例如在数值模拟或有限元素分析中。因为一旦 \(A\) 被分解，就可以快速地对任意的 \(b\) 进行求解。

在实际应用中，直接计算 \(A\) 的LU分解可能会遇到数值稳定性问题，特别是当 \(A\) 不是正定矩阵或接近奇异矩阵时。为了解决这个问题，通常会使用带部分选主元的LU分解（Pivoted LU decomposition），即在分解过程中引入一个置换矩阵 \(P\)，以保持数值稳定性，这时分解表示为 \(PA = LU\)。

R语言中的lu()函数可以用来进行LU分解，但需要注意的是，R的lu()函数返回的是一个置换矩阵 \(P\)、一个下三角矩阵 \(L\)、以及一个上三角矩阵 \(U\)，即 \(PA = LU\)。

简化的版本，用于展示基本步骤：

# 定义一个矩阵
A <- matrix(c(4, 3, -5,
              0, -2, 1,
              -1, 2, 2), byrow = TRUE, nrow = 3)

# 初始化L和U矩阵
L <- diag(1, nrow(A))  # 创建一个单位矩阵作为L的初始值
U <- A                 # 将A的初始值赋给U

# 获取矩阵的维度
n <- ncol(A)

# 执行LU分解
for (i in 1:(n-1)) {
    for (j in (i+1):n) {
        L[j, i] <- U[j, i] / U[i, i]  # 计算L矩阵的元素
        U[j, ] <- U[j, ] - L[j, i] * U[i, ]  # 更新U矩阵的行
    }
}

# 输出L和U矩阵
L

      [,1]   [,2] [,3]
[1,]  1.00  0.000    0
[2,]  0.00  1.000    0
[3,] -0.25 -1.375    1

U

     [,1] [,2]   [,3]
[1,]    4    3 -5.000
[2,]    0   -2  1.000
[3,]    0    0  2.125

# 对比结果
lu(A)

$L
      [,1]   [,2] [,3]
[1,]  1.00  0.000    0
[2,]  0.00  1.000    0
[3,] -0.25 -1.375    1

$U
     [,1] [,2]   [,3]
[1,]    4    3 -5.000
[2,]    0   -2  1.000
[3,]    0    0  2.125

在这段代码中，首先创建了单位矩阵L和一个复制自 A 的矩阵 U。对于每个元素 U[j, i]，计算它和主对角线上元素 U[i, i] 的比，这个比值就是 L 矩阵在 (j, i) 位置上的值。然后，从 U 矩阵的第 j 行中减去第 i 行乘以这个比值，以此来更新 U 矩阵。

以上代码是一个非常简化的 LU 分解过程，仅用于演示目的。它没有考虑一些特殊情况，例如主对角线上的元素为零（这会导致除以零的错误），或者矩阵 A 不是满秩矩阵（这将导致 LU 分解不存在）。在实际情况中，LU 分解可能需要进行行交换（部分主元法）来保证算法的稳定性。这样的算法实现要复杂得多，通常是由数值线性代数库提供的。在 R 中，您可以直接使用Matrix包中的lu()函数，或者pracma包中的lu()函数来得到完整的 LU 分解，包括置换矩阵。

奇异值分解

手动实现一个完整的奇异值分解（SVD）算法对我来说还办不到。

SVD 算法涉及到高级的矩阵操作，如计算矩阵的特征值和特征向量，以及进行正交化和数值优化。这些都是数值分析和矩阵理论领域中的高级主题。在计算机上，SVD 通常是通过数值算法来实现的，比如QR算法或者分解迭代方法。这些算法在实际中由专业的数值线性代数库来实现，如 LAPACK 或者 BLAS，它们经过了优化以确保计算的稳定性和效率。

尽管如此，可以整一个算法概述，这个概述是基于构造 A^TA（或者 AAT，取决于A的形状）的特征值分解，但请注意这不是实际数值软件使用的方法。

1. 对于矩阵 A（大小为m×n），计算 A^TA。
2. 计算 A^TA 的特征值和特征向量。
3. 特征向量形成了 V 矩阵，特征值的平方根形成了Σ矩阵的对角线元素。
4. 计算 U 矩阵，使得 U = AVΣ^{-1}。

在R中，计算特征值和特征向量可以使用eigen()函数，但手动执行整个SVD过程（包括确保所有数值稳定性和正确的排序）是不现实的。即使如此，下面的R代码给出了一个简化的流程：

# 假设A是一个m x n矩阵
A <- matrix(c(1, 2, 3, 4, 5, 6), nrow=3, byrow=TRUE)

# 计算A的转置乘以A
AtA <- t(A) %*% A

# 使用eigen函数计算特征值和特征向量
eig <- eigen(AtA)

# 特征向量就是V
V <- eig$vectors

# 特征值的平方根就是奇异值，构成Σ矩阵
Sigma <- sqrt(diag(eig$values))

# 计算U矩阵
U <- A %*% V %*% solve(Sigma)

# 因为特征值可能会有负值，对于奇异值分解需要保留正值
positive_singular_values <- eig$values > 0
Sigma <- sqrt(diag(eig$values[positive_singular_values]))
V <- eig$vectors[, positive_singular_values]

# 此时U的列可能不是正交的，需要进行正交化
U <- qr.Q(qr(A %*% V))

# 输出结果
U

           [,1]       [,2]
[1,] -0.2298477  0.8834610
[2,] -0.5247448  0.2407825
[3,] -0.8196419 -0.4018960

Sigma

         [,1]      [,2]
[1,] 9.525518 0.0000000
[2,] 0.000000 0.5143006

V

          [,1]       [,2]
[1,] 0.6196295 -0.7848945
[2,] 0.7848945  0.6196295

svd(A)

$d
[1] 9.5255181 0.5143006

$u
           [,1]       [,2]
[1,] -0.2298477  0.8834610
[2,] -0.5247448  0.2407825
[3,] -0.8196419 -0.4018960

$v
           [,1]       [,2]
[1,] -0.6196295 -0.7848945
[2,] -0.7848945  0.6196295

计算特征值和特征向量

手动实现计算特征值和特征向量的算法在 R 中是可能的，但它是一个非常高级的任务，通常会涉及到数值分析和矩阵理论的深入知识。计算特征值和特征向量的常用算法包括幂方法（Power Method）、反幂方法（Inverse Power Method）、QR 算法等。这些方法涉及到迭代计算，直到收敛到足够精确的值。

幂方法

使用幂方法作为一个示例，它是计算矩阵主要特征值的一种方法。幂方法可以找到模最大的特征值，但它不会找到所有特征值，也不会找到对应的特征向量。这里是一个幂方法的基本实现：

power_method <- function(A, maxiter = 1000, tol = 1e-10) {
  n <- nrow(A)
  b_k <- rep(1, n) # 初始化一个随机向量
  for (i in 1:maxiter) {
    # 通过矩阵A与向量b_k相乘计算下一个向量
    b_k1 <- A %*% b_k
    
    # 用下一个向量除以其范数来规范化
    b_k1 <- b_k1 / sqrt(sum(b_k1^2))
    
    # 计算当前向量和下一个向量的差异
    err <- sqrt(sum((b_k - b_k1)^2))
    
    # 如果当前向量和下一个向量之间的差异足够小，则停止迭代
    if (err < tol) {
      break
    }
    
    b_k <- b_k1
  }
  
  # 计算特征值
  lambda <- t(b_k) %*% A %*% b_k
  
  return(list(lambda = lambda, eigenvector = b_k))
}

# 使用一个示例矩阵
A <- matrix(c(2, 0, 0,
              0, 3, 4,
              0, 4, 9), nrow = 3, byrow = TRUE)

# 应用幂方法
result <- power_method(A)
result

$lambda
     [,1]
[1,]   11

$eigenvector
             [,1]
[1,] 3.215778e-11
[2,] 4.472136e-01
[3,] 8.944272e-01

eigen(A)

eigen() decomposition
$values
[1] 11  2  1

$vectors
          [,1] [,2]       [,3]
[1,] 0.0000000    1  0.0000000
[2,] 0.4472136    0  0.8944272
[3,] 0.8944272    0 -0.4472136

这个简单的幂方法实现可以估计矩阵 A 的最大特征值及其对应的特征向量。maxiter 参数控制迭代的最大次数，而 tol 参数定义了收敛的容差。这种方法可能不会收敛到正确的特征值，如果矩阵有多个具有相同模的特征值，或者迭代次数不够。

对于更复杂的矩阵，可能需要使用更高级的算法，如 QR 算法。幂方法仅仅能够估计矩阵的主要特征值及其对应的特征向量。对于完整的特征值和特征向量集，我们需要使用更复杂的算法，如 QR 算法或者雅可比旋转法（Jacobi’s method）等。

QR 算法

简化版本的 QR 算法。QR算法通过不断迭代应用QR分解逐步逼近矩阵的特征值。这里有一个示例：

qr_algorithm <- function(A, maxiter = 100, tol = 1e-10) {
  n <- ncol(A)
  Ak <- A
  Qk <- diag(n)
  
  for (i in 1:maxiter) {
    # QR分解
    Q <- qr.Q(qr(Ak))
    R <- qr.R(qr(Ak))
    
    # 更新A为下一个迭代的值
    Ak <- R %*% Q
    
    # 累积Q以估计特征向量
    Qk <- Qk %*% Q
    
    # 检查收敛性，这里使用的是所有非对角线元素的和
    off_diagonal_sum <- sum(abs(Ak[row(Ak) != col(Ak)]))
    if (off_diagonal_sum < tol) {
      break
    }
  }
  
  # 对角线元素是特征值
  eigenvalues <- diag(Ak)
  
  # 特征向量是累积的Q矩阵的列
  eigenvectors <- Qk
  
  return(list(values = eigenvalues, vectors = eigenvectors))
}

# 使用一个示例矩阵
A <- matrix(c(12, -51, 4,
               6, 167, -68,
               -4, 24, -41), byrow = TRUE, nrow = 3)

# 应用QR算法
result <- qr_algorithm(A)
result$values # 打印特征值

[1] 156.13668 -34.19668  16.05999

result$vectors # 打印特征向量

           [,1]       [,2]       [,3]
[1,] -0.3281474 0.37540398  0.8668282
[2,]  0.9368814 0.01207139  0.3494390
[3,]  0.1207170 0.92678268 -0.3556702

eigen(A)

eigen() decomposition
$values
[1] 156.13668 -34.19668  16.05999

$vectors
           [,1]      [,2]       [,3]
[1,]  0.3281474 0.2547577 -0.9905263
[2,] -0.9368814 0.3027934  0.0871754
[3,] -0.1207170 0.9183761  0.1061044

这段代码会给出一个矩阵所有特征值的近似值，但是不会计算特征向量。为了得到特征向量，需要进一步迭代和提取操作，这通常是由数值线性代数库如 LAPACK 来完成的。

QR 分解的原理

QR 分解的原理基于正交化过程，主要是通过 Gram-Schmidt 正交化方法或者其它更稳定的数值方法如 Householder 变换或 Givens 旋转。以下是 Gram-Schmidt 方法的基本概念：

1. Gram-Schmidt正交化:
   • 假设你有一个矩阵 \(A\) 由列向量 \(\{a_1, a_2, ..., a_n\}\) 组成。
   • 你从第一个向量 \(a_1\) 开始，将其归一化成 \(q_1\)，这是正交矩阵 \(Q\) 的第一列。
   • 接下来，从 \(a_2\) 减去它在 \(q_1\) 上的投影，得到 \(a_2\) 与 \(q_1\) 正交的向量，再归一化为 \(q_2\)。
   • 对所有的 \(a_i\) 重复这个过程，每次都从 \(a_i\) 中减去它在所有前面的 \(q_j\) 上的投影，然后归一化。
   • 这样，你就得到了一组彼此正交的单位向量 \(\{q_1, q_2, ..., q_n\}\)，它们构成 \(Q\) 的列。
2. 构造上三角矩阵 \(R\):
   • 在正交化的同时，你可以记录每一步中减去的投影的大小，这些大小构成上三角矩阵 \(R\) 的元素。
   • 具体来说，\(R\) 的对角线元素 \(r_{ii}\) 是原矩阵 \(A\) 中的 \(a_i\) 向量长度，在正交化过程中变成 \(q_i\) 向量的长度。
   • 上三角矩阵 \(R\) 的非对角线元素 \(r_{ij}\) (其中 \(i < j\)) 是向量 \(a_j\) 在向量 \(q_i\) 上的投影的长度。
3. Householder变换:
   • Householder 变换是一种更稳定的正交化方法。它通过构造反射矩阵来一次性地将一个向量的所有下方分量变为0，而不是逐个正交化。
   • 这种方法可以减少因为数值计算引起的不稳定性，并且通常用于实际的数值软件中。
4. Givens旋转:
   • Givens 旋转是另一种用于QR分解的方法，它通过旋转来将矩阵的特定元素变为0。Givens旋转特别适用于稀疏矩阵。

在实际数值计算中，Householder 变换和 Givens 旋转比 Gram-Schmidt 更为常用，因为它们在面对舍入误差时更加稳定。在计算机实现中，Householder 变换是最常用的方法，因为它为大多数情况提供了好的速度和稳定性平衡。

工作环境

devtools::session_info()

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.1 (2023-06-16 ucrt)
 os       Windows 11 x64 (build 22621)
 system   x86_64, mingw32
 ui       RTerm
 language (EN)
 collate  Chinese (Simplified)_China.utf8
 ctype    Chinese (Simplified)_China.utf8
 tz       Asia/Hong_Kong
 date     2023-11-10
 pandoc   3.1.9 @ C:/Users/HANWAN~1/AppData/Local/Pandoc/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 cachem        1.0.8   2023-05-01 [1] CRAN (R 4.3.1)
 callr         3.7.3   2022-11-02 [1] CRAN (R 4.3.1)
 cli           3.6.1   2023-03-23 [1] CRAN (R 4.3.1)
 crayon        1.5.2   2022-09-29 [1] CRAN (R 4.3.1)
 devtools      2.4.5   2022-10-11 [1] CRAN (R 4.3.2)
 digest        0.6.33  2023-07-07 [1] CRAN (R 4.3.1)
 ellipsis      0.3.2   2021-04-29 [1] CRAN (R 4.3.1)
 evaluate      0.21    2023-05-05 [1] CRAN (R 4.3.1)
 fastmap       1.1.1   2023-02-24 [1] CRAN (R 4.3.1)
 fs            1.6.3   2023-07-20 [1] CRAN (R 4.3.1)
 glue          1.6.2   2022-02-24 [1] CRAN (R 4.3.1)
 htmltools     0.5.5   2023-03-23 [1] CRAN (R 4.3.1)
 htmlwidgets   1.6.2   2023-03-17 [1] CRAN (R 4.3.1)
 httpuv        1.6.11  2023-05-11 [1] CRAN (R 4.3.1)
 jsonlite      1.8.7   2023-06-29 [1] CRAN (R 4.3.1)
 knitr         1.43    2023-05-25 [1] CRAN (R 4.3.1)
 later         1.3.1   2023-05-02 [1] CRAN (R 4.3.1)
 lifecycle     1.0.3   2022-10-07 [1] CRAN (R 4.3.1)
 magrittr      2.0.3   2022-03-30 [1] CRAN (R 4.3.1)
 memoise       2.0.1   2021-11-26 [1] CRAN (R 4.3.1)
 mime          0.12    2021-09-28 [1] CRAN (R 4.3.0)
 miniUI        0.1.1.1 2018-05-18 [1] CRAN (R 4.3.1)
 pkgbuild      1.4.2   2023-06-26 [1] CRAN (R 4.3.1)
 pkgload       1.3.2.1 2023-07-08 [1] CRAN (R 4.3.1)
 pracma      * 2.4.4   2023-11-10 [1] CRAN (R 4.3.1)
 prettyunits   1.1.1   2020-01-24 [1] CRAN (R 4.3.1)
 processx      3.8.2   2023-06-30 [1] CRAN (R 4.3.1)
 profvis       0.3.8   2023-05-02 [1] CRAN (R 4.3.1)
 promises      1.2.0.1 2021-02-11 [1] CRAN (R 4.3.1)
 ps            1.7.5   2023-04-18 [1] CRAN (R 4.3.1)
 purrr         1.0.2   2023-08-10 [1] CRAN (R 4.3.2)
 R6            2.5.1   2021-08-19 [1] CRAN (R 4.3.1)
 Rcpp          1.0.11  2023-07-06 [1] CRAN (R 4.3.1)
 remotes       2.4.2.1 2023-07-18 [1] CRAN (R 4.3.1)
 rlang         1.1.1   2023-04-28 [1] CRAN (R 4.3.1)
 rmarkdown     2.23    2023-07-01 [1] CRAN (R 4.3.1)
 rstudioapi    0.15.0  2023-07-07 [1] CRAN (R 4.3.1)
 sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.3.1)
 shiny         1.7.4.1 2023-07-06 [1] CRAN (R 4.3.1)
 stringi       1.7.12  2023-01-11 [1] CRAN (R 4.3.0)
 stringr       1.5.0   2022-12-02 [1] CRAN (R 4.3.1)
 urlchecker    1.0.1   2021-11-30 [1] CRAN (R 4.3.1)
 usethis       2.2.2   2023-07-06 [1] CRAN (R 4.3.1)
 vctrs         0.6.3   2023-06-14 [1] CRAN (R 4.3.1)
 xfun          0.39    2023-04-20 [1] CRAN (R 4.3.1)
 xtable        1.8-4   2019-04-21 [1] CRAN (R 4.3.1)
 yaml          2.3.7   2023-01-23 [1] CRAN (R 4.3.0)

 [1] C:/Users/Han Wang/AppData/Local/R/win-library/4.3
 [2] C:/Program Files/R/R-4.3.1/library

──────────────────────────────────────────────────────────────────────────────