LU矩阵分解

什么是LU分解？

LU分解是指对一个方阵 $𝐴$ 进行因式分解，通过适当的行和/或列的重排或置换，将其分解为两个因子：

• 一个下三角矩阵 $𝐿$
• 一个上三角矩阵 $𝑈$

$𝐴 = 𝐿𝑈$

1A <- matrix(c(1, -1, 2, 2), 
2            nrow = 2, 
3            byrow = TRUE)
4
5fA <- LU(A = A) # 计算矩阵 A 的 LU 分解
6A_ <- fA$P %*% fA$L %*% fA$U # 通过相乘所有因子恢复原始矩阵 A
7
8fA
9A_

说明

分解得到的 LU 因子如下：

1fA <- LU(A = A)

---

P、L、U 因子 $$\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \quad\quad L= \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix} \quad\quad U= \begin{bmatrix} 1 & -1 \\ 0 & 4 \end{bmatrix}$$

在 R 中，LU 分解还会返回置换矩阵 $P$ 以及 $L$ 和 $U$ 矩阵。为了验证分解是否正确，可以通过将所有分解因子（包括 $P$）相乘来恢复原始方阵 $A$。

1A_ <- fA$P %*% fA$L %*% fA$U

原始矩阵 A $$\begin{bmatrix} 1 & -1 \\ 2 & 2 \end{bmatrix}$$

LU分解算法

下面的 Rcpp 代码实现了 Doolittle 算法。该算法提供了一种不使用高斯消元法的替代方法，用于将矩阵 $A$ 分解为 LU 形式。该算法系统地求解 $L$ 和 $U$ 矩阵的元素，基于 $A = LU$ 所需的乘法运算结果。对于 $L$ 矩阵的元素，使用以下公式：

$$L_{ik} = A_{ik} - \displaystyle\sum_{j=0}^{k-1} L_{ij} U_{jk}$$

对于 $U$ 矩阵的元素，使用以下公式：

$$U_{ik} = \cfrac{A_{ik} - \displaystyle\sum_{j=0}^{k-1} L_{ij} U_{jk}}{L_{kk}}$$

上述公式可以在 R 中如下实现：

1Rcpp::sourceCpp('/usercode/main.cpp')
2
3A <- matrix(c(1, -1, 2, 2), 
4            nrow = 2, 
5            byrow = TRUE)
6
7decoA <- LU_fact(A) # 执行 LU 分解
8A_retrieved <- decoA$L_mat %*% decoA$U_mat # 通过相乘 L 因子和 U 因子恢复 A
9all.equal(A_retrieved, A)

在 main.r 文件中：

• 第 7 行：使用 R 中的 LU_fact() 函数对矩阵 A 进行 LU 分解。
• 第 8 行：通过将 L 因子 L_mat 与 U 因子 U_mat 相乘，恢复矩阵 $A$。
• 第 9 行：将恢复出的矩阵 A_retrieved 与原始矩阵 A 进行比较。
• 最后，将原始方阵与 L 矩阵和 U 矩阵的乘积进行对比，以确认矩阵的 LU 分解是正确的。

1decoA <- LU_fact(A)
2A_retrieved <- decoA$L_mat %*% decoA$U_mat 
3all.equal(A_retrieved, A)
4
5## [1] TRUE

使用 Armadillo 的 LU 分解

Armadillo 库提供了 lu() 函数用于 LU 分解。$L$ 和 $U$ 矩阵以引用方式传入 lu() 函数。

1Rcpp::sourceCpp('/usercode/main.cpp')
2
3A <- matrix(c(1, -1, 2, 2), 
4            nrow = 2, 
5            byrow = TRUE)
6
7decoA_A <- LU_fact_A(A)
8A_retrieved_arma <- decoA_A$L_mat %*% decoA_A$U_mat 
9all.equal(A_retrieved_arma, A)

将原始方阵 $A$ 与使用 Armadillo 函数计算出的 $L$ 矩阵和 $U$ 矩阵的乘积进行比较，确认 LU 分解正确。

1decoA_A <- LU_fact_A(A)
2A_retrieved_arma <- decoA_A$L_mat %*% decoA_A$U_mat 
3all.equal(A_retrieved_arma, A)
4
5## [1] TRUE

使用 Eigen 的 LU 分解

在 Eigen 中获取 LU 分解因子稍显复杂，需要引入以下类：

• 使用 PartialPivLU 类进行分解计算；
• 使用 UpLoType 类提取下三角和上三角部分（此处下三角为单位下三角）。

1Rcpp::sourceCpp('/usercode/main.cpp')
2
3A <- matrix(c(1, -1, 2, 2), 
4            nrow = 2, 
5            byrow = TRUE)
6
7decoA_E <- LU_fact_E(A) # 执行 LU 分解
8A_retrieved_Eigen <- decoA_E$perm %*% decoA_E$L_mat %*% decoA_E$U_mat # 通过相乘 L 和 U 因子恢复 A
9all.equal(A_retrieved_Eigen, A)

将原始方阵 $A$ 与使用 Eigen 函数计算出的 $L$ 矩阵和 $U$ 矩阵的乘积进行比较，确认分解正确。

1decoA_E <- LU_fact_E(A)
2A_retrieved_Eigen <- decoA_E$perm %*% decoA_E$L_mat %*% decoA_E$U_mat 
3all.equal(A_retrieved_Eigen, A)
4
5## [1] TRUE