编写一个可跟踪的 R 函数，模拟 LAPACK 的 dgetrf 以进行 LU 分解
Write a trackable R function that mimics LAPACK's dgetrf for LU factorization
R 核心中没有 LU 分解函数。尽管这种因式分解是 solve 的一个步骤，但它并未作为独立函数明确提供。我们可以为此编写一个 R 函数吗？它需要模仿 LAPACK 例程 dgetrf. Matrix package has an lu function 这很好，但如果我们可以编写一个 trackable R 函数会更好，它可以
- 分解矩阵直到某一列/行和return中间结果；
- 继续从中间结果分解到另一列/行或最后。
此功能对于教育和调试目的都非常有用。教育的好处是显而易见的，因为我们可以逐列说明因式分解/高斯消元法。为了调试使用，这里有两个例子。
在中询问为什么R中的LU因式分解和Python给出不同的结果。我们可以清楚地看到两个软件 return 相同的第一个枢轴和第二个枢轴，但不是第三个枢轴。所以当分解进行到第 3 行/列时一定有一些有趣的事情。如果我们能够检索到该临时结果以进行调查，那就太好了。
在 中，LU 分解对于这种类型的矩阵是不稳定的。在我的回答中，给出了一个 3 x 3 矩阵作为示例。我希望 solve 产生错误抱怨 U[3, 3] = 0，但是 运行 solve 有几次我发现 solve 有时会成功。因此，对于数值调查，我想知道当分解进行到第二列/行时会发生什么。
由于该函数是用纯R代码编写的，对于中等到大的矩阵，预计速度会很慢。但是性能不是问题，至于教育和调试我们只用了一个小矩阵。
dgetrf 简介
LAPACK 的 dgetrf 使用行旋转计算 LU 分解：A = PLU。在因式分解退出时，
L为单位下三角矩阵，存放在A的下三角部分；U为上三角矩阵，存储在A的上三角部分；P 是一个行置换矩阵，存储为单独的置换索引向量。
除非主元正好为零（不符合某些公差），否则应继续分解。
我从什么开始
编写既不使用行旋转也不使用 "pause / continue" 选项的 LU 分解并不具有挑战性：
LU <- function (A) {
## check dimension
n <- dim(A)
if (n[1] != n[2]) stop("'A' must be a square matrix")
n <- n[1]
## Gaussian elimination
for (j in 1:(n - 1)) {
ind <- (j + 1):n
## check if the pivot is EXACTLY 0
piv <- A[j, j]
if (piv == 0) stop(sprintf("system is exactly singular: U[%d, %d] = 0", j, j))
l <- A[ind, j] / piv
## update `L` factor
A[ind, j] <- l
## update `U` factor by Gaussian elimination
A[ind, ind] <- A[ind, ind] - tcrossprod(l, A[j, ind])
}
A
}
这显示在不需要旋转时给出正确的结果：
A <- structure(c(0.923065107548609, 0.922819485189393, 0.277002309216186,
0.532856695353985, 0.481061384081841, 0.0952619954477996,
0.261916425777599, 0.433514681644738, 0.677919807843864,
0.771985625848174, 0.705952850636095, 0.873727774480358,
0.28782021952793, 0.863347264472395, 0.627262107795104,
0.187472499441355), .Dim = c(4L, 4L))
oo <- LU(A)
oo
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.9997339 -0.3856714 0.09424621 0.5756036
#[3,] 0.3000897 -0.3048058 0.53124291 0.7163376
#[4,] 0.5772688 -0.4040044 0.97970570 -0.4479307
L <- diag(4)
low <- lower.tri(L)
L[low] <- oo[low]
L
# [,1] [,2] [,3] [,4]
#[1,] 1.0000000 0.0000000 0.0000000 0
#[2,] 0.9997339 1.0000000 0.0000000 0
#[3,] 0.3000897 -0.3048058 1.0000000 0
#[4,] 0.5772688 -0.4040044 0.9797057 1
U <- oo
U[low] <- 0
U
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.0000000 -0.3856714 0.09424621 0.5756036
#[3,] 0.0000000 0.0000000 0.53124291 0.7163376
#[4,] 0.0000000 0.0000000 0.00000000 -0.4479307
与 Matrix 包中的 lu 比较：
library(Matrix)
rr <- expand(lu(A))
rr
#$L
#4 x 4 Matrix of class "dtrMatrix" (unitriangular)
# [,1] [,2] [,3] [,4]
#[1,] 1.0000000 . . .
#[2,] 0.9997339 1.0000000 . .
#[3,] 0.3000897 -0.3048058 1.0000000 .
#[4,] 0.5772688 -0.4040044 0.9797057 1.0000000
#
#$U
#4 x 4 Matrix of class "dtrMatrix"
# [,1] [,2] [,3] [,4]
#[1,] 0.92306511 0.48106138 0.67791981 0.28782022
#[2,] . -0.38567138 0.09424621 0.57560363
#[3,] . . 0.53124291 0.71633755
#[4,] . . . -0.44793070
#
#$P
#4 x 4 sparse Matrix of class "pMatrix"
#
#[1,] | . . .
#[2,] . | . .
#[3,] . . | .
#[4,] . . . |
现在考虑排列 A:
B <- A[c(4, 3, 1, 2), ]
LU(B)
# [,1] [,2] [,3] [,4]
#[1,] 0.5328567 0.43351468 0.8737278 0.1874725
#[2,] 0.5198439 0.03655646 0.2517508 0.5298057
#[3,] 1.7322952 -7.38348421 1.0231633 3.8748743
#[4,] 1.7318343 -17.93154011 3.6876940 -4.2504433
结果不同于LU(A)。但是，由于 Matrix::lu 执行行旋转，因此 lu(B) 的结果仅在置换矩阵中与 lu(A) 不同：
expand(lu(B))$P
#4 x 4 sparse Matrix of class "pMatrix"
#
#[1,] . . . |
#[2,] . . | .
#[3,] | . . .
#[4,] . | . .
让我们一一添加这些功能。
行旋转
这并不太难。
假设 A 是 n x n。初始化一个排列索引向量pivot <- 1:n。在第 j 列，我们扫描 A[j:n, j] 以获得最大绝对值。假设是A[m, j]。如果 m > j 我们进行行交换 A[m, ] <-> A[j, ]。同时我们做一个排列pivot[j] <-> pivot[m]。主元后消元法与不主元分解相同，所以我们可以重用LU.
函数的代码
LUP <- function (A) {
## check dimension
n <- dim(A)
if (n[1] != n[2]) stop("'A' must be a square matrix")
n <- n[1]
## LU factorization from the beginning to the end
from <- 1
to <- (n - 1)
pivot <- 1:n
## Gaussian elimination
for (j in from:to) {
## select pivot
m <- which.max(abs(A[j:n, j]))
## A[j - 1 + m, j] is the pivot
if (m > 1L) {
## row exchange
tmp <- A[j, ]; A[j, ] <- A[j - 1 + m, ]; A[j - 1 + m, ] <- tmp
tmp <- pivot[j]; pivot[j] <- pivot[j - 1 + m]; pivot[j - 1 + m] <- tmp
}
ind <- (j + 1):n
## check if the pivot is EXACTLY 0
piv <- A[j, j]
if (piv == 0) {
stop(sprintf("system is exactly singular: U[%d, %d] = 0", j, j))
}
l <- A[ind, j] / piv
## update `L` factor
A[ind, j] <- l
## update `U` factor by Gaussian elimination
A[ind, ind] <- A[ind, ind] - tcrossprod(l, A[j, ind])
}
## add `pivot` as an attribute and return `A`
structure(A, pivot = pivot)
}
在问题中尝试矩阵 B，LUP(B) 与 LU(A) 相同，但有一个额外的置换索引向量。
oo <- LUP(B)
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.9997339 -0.3856714 0.09424621 0.5756036
#[3,] 0.3000897 -0.3048058 0.53124291 0.7163376
#[4,] 0.5772688 -0.4040044 0.97970570 -0.4479307
#attr(,"pivot")
#[1] 3 4 2 1
这里是提取L,U,P:
的效用函数
exLUP <- function (LUPftr) {
L <- diag(1, nrow(LUPftr), ncol(LUPftr))
low <- lower.tri(L)
L[low] <- LUPftr[low]
U <- LUPftr[1:nrow(LUPftr), ] ## use "[" to drop attributes
U[low] <- 0
list(L = L, U = U, P = attr(LUPftr, "pivot"))
}
rr <- exLUP(oo)
#$L
# [,1] [,2] [,3] [,4]
#[1,] 1.0000000 0.0000000 0.0000000 0
#[2,] 0.9997339 1.0000000 0.0000000 0
#[3,] 0.3000897 -0.3048058 1.0000000 0
#[4,] 0.5772688 -0.4040044 0.9797057 1
#
#$U
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.0000000 -0.3856714 0.09424621 0.5756036
#[3,] 0.0000000 0.0000000 0.53124291 0.7163376
#[4,] 0.0000000 0.0000000 0.00000000 -0.4479307
#
#$P
#[1] 3 4 2 1
注意排列索引returned真的是为了PA = LU（可能是教科书上用得最多的）：
all.equal( B[rr$P, ], with(rr, L %*% U) )
#[1] TRUE
要获得由 LAPACK 编辑的 return 排列索引，即 A = PLU 中的排列索引，请执行 order(rr$P).
all.equal( B, with(rr, (L %*% U)[order(P), ]) )
#[1] TRUE
“暂停/继续”选项
添加“暂停/继续”功能有点棘手，因为我们需要一些方法来记录不完全分解停止的位置，以便我们稍后可以从那里获取它。
假设我们要将功能 LUP 增强为一个新功能 LUP2。考虑添加一个参数 to。因式分解将在 A[to, to] 完成后停止，并将与 A[to + 1, to + 1] 一起使用。我们可以将此 to 以及临时 pivot 向量存储为 A 和 return 的属性。稍后当我们将这个临时结果传回LUP2时，需要先检查这些属性是否存在。如果是这样，它知道应该从哪里开始；否则它只是从头开始。
LUP2 <- function (A, to = NULL) {
## check dimension
n <- dim(A)
if (n[1] != n[2]) stop("'A' must be a square matrix")
n <- n[1]
## ensure that "to" has a valid value
## if it is not provided, set it to (n - 1) so that we complete factorization of `A`
## if provided, it can not be larger than (n - 1); otherwise it is reset to (n - 1)
if (is.null(to)) to <- n - 1L
else if (to > n - 1L) {
warning(sprintf("provided 'to' too big; reset to maximum possible value: %d", n - 1L))
to <- n - 1L
}
## is `A` an intermediate result of a previous, unfinished LU factorization?
## if YES, it should have a "to" attribute, telling where the previous factorization stopped
## if NO, a new factorization starting from `A[1, 1]` is performed
from <- attr(A, "to")
if (!is.null(from)) {
## so we continue factorization, but need to make sure there is work to do
from <- from + 1L
if (from >= n) {
warning("LU factorization of is already completed; return input as it is")
return(A)
}
if (from > to) {
stop(sprintf("please provide a bigger 'to' between %d and %d", from, n - 1L))
}
## extract "pivot"
pivot <- attr(A, "pivot")
} else {
## we start a new factorization
from <- 1
pivot <- 1:n
}
## LU factorization from `A[from, from]` to `A[to, to]`
## the following code reuses function `LUP`'s code
for (j in from:to) {
## select pivot
m <- which.max(abs(A[j:n, j]))
## A[j - 1 + m, j] is the pivot
if (m > 1L) {
## row exchange
tmp <- A[j, ]; A[j, ] <- A[j - 1 + m, ]; A[j - 1 + m, ] <- tmp
tmp <- pivot[j]; pivot[j] <- pivot[j - 1 + m]; pivot[j - 1 + m] <- tmp
}
ind <- (j + 1):n
## check if the pivot is EXACTLY 0
piv <- A[j, j]
if (piv == 0) {
stop(sprintf("system is exactly singular: U[%d, %d] = 0", j, j))
}
l <- A[ind, j] / piv
## update `L` factor
A[ind, j] <- l
## update `U` factor by Gaussian elimination
A[ind, ind] <- A[ind, ind] - tcrossprod(l, A[j, ind])
}
## update attributes of `A` and return `A`
structure(A, to = to, pivot = pivot)
}
尝试在问题中使用矩阵 B。假设我们想在处理完 2 列/行后停止分解。
oo <- LUP2(B, 2)
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.9997339 -0.3856714 0.09424621 0.5756036
#[3,] 0.5772688 -0.4040044 0.52046170 0.2538693
#[4,] 0.3000897 -0.3048058 0.53124291 0.7163376
#attr(,"to")
#[1] 2
#attr(,"pivot")
#[1] 3 4 1 2
由于分解不完全，U 因子不是上三角。这是提取它的辅助函数。
## usable for all functions: `LU`, `LUP` and `LUP2`
## for `LUP2` the attribute "to" is used;
## for other two we can simply zero the lower triangular of `A`
getU <- function (A) {
attr(A, "pivot") <- NULL
to <- attr(A, "to")
if (is.null(to)) {
A[lower.tri(A)] <- 0
} else {
n <- nrow(A)
len <- (n - 1):(n - to)
zero_ind <- sequence(len)
offset <- seq.int(1L, by = n + 1L, length = to)
zero_ind <- zero_ind + rep.int(offset, len)
A[zero_ind] <- 0
}
A
}
getU(oo)
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.0000000 -0.3856714 0.09424621 0.5756036
#[3,] 0.0000000 0.0000000 0.52046170 0.2538693
#[4,] 0.0000000 0.0000000 0.53124291 0.7163376
#attr(,"to")
#[1] 2
现在我们可以继续因式分解了：
LUP2(oo, 1)
#Error in LUP2(oo, 1) : please provide a bigger 'to' between 3 and 3
糟糕，我们错误地将一个不可行的值 to = 1 传递给了 LUP2，因为临时结果已经处理了 2 列/行并且无法撤消。该函数告诉我们只能向前移动，将 to 设置为 3 到 3 之间的任意整数。如果我们传入大于 3 的值，则会产生警告并将 to 重置为最大可能值。
oo <- LUP2(oo, 10)
#Warning message:
#In LUP2(oo, 10) :
# provided 'to' too big; reset to maximum possible value: 3
我们有 U 因素
getU(oo)
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.0000000 -0.3856714 0.09424621 0.5756036
#[3,] 0.0000000 0.0000000 0.53124291 0.7163376
#[4,] 0.0000000 0.0000000 0.00000000 -0.4479307
#attr(,"to")
#[1] 3
oo 现在是完全分解结果。如果我们仍然要求 LUP2 更新呢？
## without providing "to", it defaults to factorize till the end
oo <- LUP2(oo)
#Warning message:
#In LUP2(oo) :
# LU factorization is already completed; return input as it is
它告诉您无法再做任何事情，并且 return 按原样输入。
最后让我们试试奇异方阵。
## this 4 x 4 matrix has rank 1
S <- tcrossprod(1:4, 2:5)
LUP2(S)
#Error in LUP2(S) : system is exactly singular: U[2, 2] = 0
## traceback
LUP2(S, to = 1)
# [,1] [,2] [,3] [,4]
#[1,] 8.00 12 16 20
#[2,] 0.50 0 0 0
#[3,] 0.75 0 0 0
#[4,] 0.25 0 0 0
#attr(,"to")
#[1] 1
#attr(,"pivot")
#[1] 4 2 3 1
R 核心中没有 LU 分解函数。尽管这种因式分解是 solve 的一个步骤，但它并未作为独立函数明确提供。我们可以为此编写一个 R 函数吗？它需要模仿 LAPACK 例程 dgetrf. Matrix package has an lu function 这很好，但如果我们可以编写一个 trackable R 函数会更好，它可以
- 分解矩阵直到某一列/行和return中间结果；
- 继续从中间结果分解到另一列/行或最后。
此功能对于教育和调试目的都非常有用。教育的好处是显而易见的，因为我们可以逐列说明因式分解/高斯消元法。为了调试使用，这里有两个例子。
在
在 solve 产生错误抱怨 U[3, 3] = 0，但是 运行 solve 有几次我发现 solve 有时会成功。因此，对于数值调查，我想知道当分解进行到第二列/行时会发生什么。
由于该函数是用纯R代码编写的，对于中等到大的矩阵，预计速度会很慢。但是性能不是问题，至于教育和调试我们只用了一个小矩阵。
dgetrf 简介
LAPACK 的 dgetrf 使用行旋转计算 LU 分解：A = PLU。在因式分解退出时，
L为单位下三角矩阵，存放在A的下三角部分；U为上三角矩阵，存储在A的上三角部分；P是一个行置换矩阵，存储为单独的置换索引向量。
除非主元正好为零（不符合某些公差），否则应继续分解。
我从什么开始
编写既不使用行旋转也不使用 "pause / continue" 选项的 LU 分解并不具有挑战性：
LU <- function (A) {
## check dimension
n <- dim(A)
if (n[1] != n[2]) stop("'A' must be a square matrix")
n <- n[1]
## Gaussian elimination
for (j in 1:(n - 1)) {
ind <- (j + 1):n
## check if the pivot is EXACTLY 0
piv <- A[j, j]
if (piv == 0) stop(sprintf("system is exactly singular: U[%d, %d] = 0", j, j))
l <- A[ind, j] / piv
## update `L` factor
A[ind, j] <- l
## update `U` factor by Gaussian elimination
A[ind, ind] <- A[ind, ind] - tcrossprod(l, A[j, ind])
}
A
}
这显示在不需要旋转时给出正确的结果：
A <- structure(c(0.923065107548609, 0.922819485189393, 0.277002309216186,
0.532856695353985, 0.481061384081841, 0.0952619954477996,
0.261916425777599, 0.433514681644738, 0.677919807843864,
0.771985625848174, 0.705952850636095, 0.873727774480358,
0.28782021952793, 0.863347264472395, 0.627262107795104,
0.187472499441355), .Dim = c(4L, 4L))
oo <- LU(A)
oo
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.9997339 -0.3856714 0.09424621 0.5756036
#[3,] 0.3000897 -0.3048058 0.53124291 0.7163376
#[4,] 0.5772688 -0.4040044 0.97970570 -0.4479307
L <- diag(4)
low <- lower.tri(L)
L[low] <- oo[low]
L
# [,1] [,2] [,3] [,4]
#[1,] 1.0000000 0.0000000 0.0000000 0
#[2,] 0.9997339 1.0000000 0.0000000 0
#[3,] 0.3000897 -0.3048058 1.0000000 0
#[4,] 0.5772688 -0.4040044 0.9797057 1
U <- oo
U[low] <- 0
U
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.0000000 -0.3856714 0.09424621 0.5756036
#[3,] 0.0000000 0.0000000 0.53124291 0.7163376
#[4,] 0.0000000 0.0000000 0.00000000 -0.4479307
与 Matrix 包中的 lu 比较：
library(Matrix)
rr <- expand(lu(A))
rr
#$L
#4 x 4 Matrix of class "dtrMatrix" (unitriangular)
# [,1] [,2] [,3] [,4]
#[1,] 1.0000000 . . .
#[2,] 0.9997339 1.0000000 . .
#[3,] 0.3000897 -0.3048058 1.0000000 .
#[4,] 0.5772688 -0.4040044 0.9797057 1.0000000
#
#$U
#4 x 4 Matrix of class "dtrMatrix"
# [,1] [,2] [,3] [,4]
#[1,] 0.92306511 0.48106138 0.67791981 0.28782022
#[2,] . -0.38567138 0.09424621 0.57560363
#[3,] . . 0.53124291 0.71633755
#[4,] . . . -0.44793070
#
#$P
#4 x 4 sparse Matrix of class "pMatrix"
#
#[1,] | . . .
#[2,] . | . .
#[3,] . . | .
#[4,] . . . |
现在考虑排列 A:
B <- A[c(4, 3, 1, 2), ]
LU(B)
# [,1] [,2] [,3] [,4]
#[1,] 0.5328567 0.43351468 0.8737278 0.1874725
#[2,] 0.5198439 0.03655646 0.2517508 0.5298057
#[3,] 1.7322952 -7.38348421 1.0231633 3.8748743
#[4,] 1.7318343 -17.93154011 3.6876940 -4.2504433
结果不同于LU(A)。但是，由于 Matrix::lu 执行行旋转，因此 lu(B) 的结果仅在置换矩阵中与 lu(A) 不同：
expand(lu(B))$P
#4 x 4 sparse Matrix of class "pMatrix"
#
#[1,] . . . |
#[2,] . . | .
#[3,] | . . .
#[4,] . | . .
让我们一一添加这些功能。
行旋转
这并不太难。
假设 A 是 n x n。初始化一个排列索引向量pivot <- 1:n。在第 j 列，我们扫描 A[j:n, j] 以获得最大绝对值。假设是A[m, j]。如果 m > j 我们进行行交换 A[m, ] <-> A[j, ]。同时我们做一个排列pivot[j] <-> pivot[m]。主元后消元法与不主元分解相同，所以我们可以重用LU.
LUP <- function (A) {
## check dimension
n <- dim(A)
if (n[1] != n[2]) stop("'A' must be a square matrix")
n <- n[1]
## LU factorization from the beginning to the end
from <- 1
to <- (n - 1)
pivot <- 1:n
## Gaussian elimination
for (j in from:to) {
## select pivot
m <- which.max(abs(A[j:n, j]))
## A[j - 1 + m, j] is the pivot
if (m > 1L) {
## row exchange
tmp <- A[j, ]; A[j, ] <- A[j - 1 + m, ]; A[j - 1 + m, ] <- tmp
tmp <- pivot[j]; pivot[j] <- pivot[j - 1 + m]; pivot[j - 1 + m] <- tmp
}
ind <- (j + 1):n
## check if the pivot is EXACTLY 0
piv <- A[j, j]
if (piv == 0) {
stop(sprintf("system is exactly singular: U[%d, %d] = 0", j, j))
}
l <- A[ind, j] / piv
## update `L` factor
A[ind, j] <- l
## update `U` factor by Gaussian elimination
A[ind, ind] <- A[ind, ind] - tcrossprod(l, A[j, ind])
}
## add `pivot` as an attribute and return `A`
structure(A, pivot = pivot)
}
在问题中尝试矩阵 B，LUP(B) 与 LU(A) 相同，但有一个额外的置换索引向量。
oo <- LUP(B)
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.9997339 -0.3856714 0.09424621 0.5756036
#[3,] 0.3000897 -0.3048058 0.53124291 0.7163376
#[4,] 0.5772688 -0.4040044 0.97970570 -0.4479307
#attr(,"pivot")
#[1] 3 4 2 1
这里是提取L,U,P:
exLUP <- function (LUPftr) {
L <- diag(1, nrow(LUPftr), ncol(LUPftr))
low <- lower.tri(L)
L[low] <- LUPftr[low]
U <- LUPftr[1:nrow(LUPftr), ] ## use "[" to drop attributes
U[low] <- 0
list(L = L, U = U, P = attr(LUPftr, "pivot"))
}
rr <- exLUP(oo)
#$L
# [,1] [,2] [,3] [,4]
#[1,] 1.0000000 0.0000000 0.0000000 0
#[2,] 0.9997339 1.0000000 0.0000000 0
#[3,] 0.3000897 -0.3048058 1.0000000 0
#[4,] 0.5772688 -0.4040044 0.9797057 1
#
#$U
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.0000000 -0.3856714 0.09424621 0.5756036
#[3,] 0.0000000 0.0000000 0.53124291 0.7163376
#[4,] 0.0000000 0.0000000 0.00000000 -0.4479307
#
#$P
#[1] 3 4 2 1
注意排列索引returned真的是为了PA = LU（可能是教科书上用得最多的）：
all.equal( B[rr$P, ], with(rr, L %*% U) )
#[1] TRUE
要获得由 LAPACK 编辑的 return 排列索引，即 A = PLU 中的排列索引，请执行 order(rr$P).
all.equal( B, with(rr, (L %*% U)[order(P), ]) )
#[1] TRUE
“暂停/继续”选项
添加“暂停/继续”功能有点棘手，因为我们需要一些方法来记录不完全分解停止的位置，以便我们稍后可以从那里获取它。
假设我们要将功能 LUP 增强为一个新功能 LUP2。考虑添加一个参数 to。因式分解将在 A[to, to] 完成后停止，并将与 A[to + 1, to + 1] 一起使用。我们可以将此 to 以及临时 pivot 向量存储为 A 和 return 的属性。稍后当我们将这个临时结果传回LUP2时，需要先检查这些属性是否存在。如果是这样，它知道应该从哪里开始；否则它只是从头开始。
LUP2 <- function (A, to = NULL) {
## check dimension
n <- dim(A)
if (n[1] != n[2]) stop("'A' must be a square matrix")
n <- n[1]
## ensure that "to" has a valid value
## if it is not provided, set it to (n - 1) so that we complete factorization of `A`
## if provided, it can not be larger than (n - 1); otherwise it is reset to (n - 1)
if (is.null(to)) to <- n - 1L
else if (to > n - 1L) {
warning(sprintf("provided 'to' too big; reset to maximum possible value: %d", n - 1L))
to <- n - 1L
}
## is `A` an intermediate result of a previous, unfinished LU factorization?
## if YES, it should have a "to" attribute, telling where the previous factorization stopped
## if NO, a new factorization starting from `A[1, 1]` is performed
from <- attr(A, "to")
if (!is.null(from)) {
## so we continue factorization, but need to make sure there is work to do
from <- from + 1L
if (from >= n) {
warning("LU factorization of is already completed; return input as it is")
return(A)
}
if (from > to) {
stop(sprintf("please provide a bigger 'to' between %d and %d", from, n - 1L))
}
## extract "pivot"
pivot <- attr(A, "pivot")
} else {
## we start a new factorization
from <- 1
pivot <- 1:n
}
## LU factorization from `A[from, from]` to `A[to, to]`
## the following code reuses function `LUP`'s code
for (j in from:to) {
## select pivot
m <- which.max(abs(A[j:n, j]))
## A[j - 1 + m, j] is the pivot
if (m > 1L) {
## row exchange
tmp <- A[j, ]; A[j, ] <- A[j - 1 + m, ]; A[j - 1 + m, ] <- tmp
tmp <- pivot[j]; pivot[j] <- pivot[j - 1 + m]; pivot[j - 1 + m] <- tmp
}
ind <- (j + 1):n
## check if the pivot is EXACTLY 0
piv <- A[j, j]
if (piv == 0) {
stop(sprintf("system is exactly singular: U[%d, %d] = 0", j, j))
}
l <- A[ind, j] / piv
## update `L` factor
A[ind, j] <- l
## update `U` factor by Gaussian elimination
A[ind, ind] <- A[ind, ind] - tcrossprod(l, A[j, ind])
}
## update attributes of `A` and return `A`
structure(A, to = to, pivot = pivot)
}
尝试在问题中使用矩阵 B。假设我们想在处理完 2 列/行后停止分解。
oo <- LUP2(B, 2)
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.9997339 -0.3856714 0.09424621 0.5756036
#[3,] 0.5772688 -0.4040044 0.52046170 0.2538693
#[4,] 0.3000897 -0.3048058 0.53124291 0.7163376
#attr(,"to")
#[1] 2
#attr(,"pivot")
#[1] 3 4 1 2
由于分解不完全，U 因子不是上三角。这是提取它的辅助函数。
## usable for all functions: `LU`, `LUP` and `LUP2`
## for `LUP2` the attribute "to" is used;
## for other two we can simply zero the lower triangular of `A`
getU <- function (A) {
attr(A, "pivot") <- NULL
to <- attr(A, "to")
if (is.null(to)) {
A[lower.tri(A)] <- 0
} else {
n <- nrow(A)
len <- (n - 1):(n - to)
zero_ind <- sequence(len)
offset <- seq.int(1L, by = n + 1L, length = to)
zero_ind <- zero_ind + rep.int(offset, len)
A[zero_ind] <- 0
}
A
}
getU(oo)
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.0000000 -0.3856714 0.09424621 0.5756036
#[3,] 0.0000000 0.0000000 0.52046170 0.2538693
#[4,] 0.0000000 0.0000000 0.53124291 0.7163376
#attr(,"to")
#[1] 2
现在我们可以继续因式分解了：
LUP2(oo, 1)
#Error in LUP2(oo, 1) : please provide a bigger 'to' between 3 and 3
糟糕，我们错误地将一个不可行的值 to = 1 传递给了 LUP2，因为临时结果已经处理了 2 列/行并且无法撤消。该函数告诉我们只能向前移动，将 to 设置为 3 到 3 之间的任意整数。如果我们传入大于 3 的值，则会产生警告并将 to 重置为最大可能值。
oo <- LUP2(oo, 10)
#Warning message:
#In LUP2(oo, 10) :
# provided 'to' too big; reset to maximum possible value: 3
我们有 U 因素
getU(oo)
# [,1] [,2] [,3] [,4]
#[1,] 0.9230651 0.4810614 0.67791981 0.2878202
#[2,] 0.0000000 -0.3856714 0.09424621 0.5756036
#[3,] 0.0000000 0.0000000 0.53124291 0.7163376
#[4,] 0.0000000 0.0000000 0.00000000 -0.4479307
#attr(,"to")
#[1] 3
oo 现在是完全分解结果。如果我们仍然要求 LUP2 更新呢？
## without providing "to", it defaults to factorize till the end
oo <- LUP2(oo)
#Warning message:
#In LUP2(oo) :
# LU factorization is already completed; return input as it is
它告诉您无法再做任何事情，并且 return 按原样输入。
最后让我们试试奇异方阵。
## this 4 x 4 matrix has rank 1
S <- tcrossprod(1:4, 2:5)
LUP2(S)
#Error in LUP2(S) : system is exactly singular: U[2, 2] = 0
## traceback
LUP2(S, to = 1)
# [,1] [,2] [,3] [,4]
#[1,] 8.00 12 16 20
#[2,] 0.50 0 0 0
#[3,] 0.75 0 0 0
#[4,] 0.25 0 0 0
#attr(,"to")
#[1] 1
#attr(,"pivot")
#[1] 4 2 3 1