LU分解中的行主元素分解

Question

下面的函数在进行LU分解时没有使用行交换。在R中是否存在一个使用行交换的LU分解函数？

> require(Matrix)
> expand(lu(matrix(rnorm(16),4,4)))
$L
4 x 4 Matrix of class "dtrMatrix"
     [,1]        [,2]        [,3]        [,4]       
[1,]  1.00000000           .           .           .
[2,]  0.13812836  1.00000000           .           .
[3,]  0.27704442  0.39877260  1.00000000           .
[4,] -0.08512341 -0.24699820  0.04347201  1.00000000

$U
4 x 4 Matrix of class "dtrMatrix"
     [,1]       [,2]       [,3]       [,4]      
[1,]  1.5759031 -0.2074224 -1.5334082 -0.5959756
[2,]          . -1.3096874 -0.6301727  1.1953838
[3,]          .          .  1.6316292  0.6256619
[4,]          .          .          .  0.8078140

$P
4 x 4 sparse Matrix of class "pMatrix"

[1,] | . . .
[2,] . | . .
[3,] . . . |
[4,] . . | .

- user236215

---

分享一个关于LU分解的教育性问答：编写一个可追踪的R函数，模仿LAPACK的dgetrf进行LU分解。 - Zheyuan Li

Answer 1

在R语言中，lu函数使用部分(行)主元消去法进行计算。由于您的例子没有给出原始矩阵，因此我将创建一个新的示例来进行演示。

在R语言中，lu函数计算A=PLU，这相当于通过排列矩阵P^-1的行来计算矩阵A的LU分解：P^-1A=LU。有关更多信息，请参见Matrix包文档。

示例：

> A <- matrix(c(1, 1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, 1, -1, 1, -1), 4)
> A
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1   -1   -1
[3,]    1   -1   -1    1
[4,]    1   -1    1   -1

这里是 "L" 因素：

> luDec <- lu(A)
> L <- expand(luDec)$L
> L
4 x 4 Matrix of class "dtrMatrix" (unitriangular)
     [,1] [,2] [,3] [,4]
[1,]    1    .    .    .
[2,]    1    1    .    .
[3,]    1    0    1    .
[4,]    1    1   -1    1

这里是“U”系数：

> U <- expand(luDec)$U
> U
4 x 4 Matrix of class "dtrMatrix"
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    .   -2   -2    0
[3,]    .    .   -2   -2
[4,]    .    .    .   -4

这是排列矩阵：

> P <- expand(luDec)$P
> P
4 x 4 sparse Matrix of class "pMatrix"

[1,] | . . .
[2,] . . | .
[3,] . | . .
[4,] . . . |

我们可以看到LU是A的行置换版本:

> L %*% U
4 x 4 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1   -1   -1    1
[3,]    1    1   -1   -1
[4,]    1   -1    1   -1

回到原始身份A = PLU，我们可以恢复A（与上面的A进行比较）：

> P %*% L %*% U
4 x 4 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1   -1   -1
[3,]    1   -1   -1    1
[4,]    1   -1    1   -1

- David Alber

---

+1 David。非常好的答案。我在想我们如何解释置换矩阵P。各个位置上的竖线和点是什么意思？ - ColorStatistics

Answer 2

或许 this 可以胜任这项工作。但是，没有 Windows 二进制文件，我无法尝试它。

- gd047

---

Magma需要底层的magma库以及一张NVidia显卡进行GPU计算。在R中有枢轴解决方案 - 我会在r-help上询问。 - Dirk Eddelbuettel