r - 编写一个可跟踪的R函数，该函数模仿LAPACK的dgetrf进行LU分解

R核心中没有LU分解功能。尽管这种分解是solve的一个步骤，但并未明确将其作为独立函数使用。我们可以为此编写一个R函数吗？它需要模拟LAPACK例程 dgetrf 。 Matrix包中有一个 lu function很好，但是如果我们可以编写一个可跟踪的 R函数会更好，它可以

将矩阵分解为一定的列/行，并返回中间结果；

从中间结果继续分解到另一列/行或结尾。

此功能对于教育和调试目的都是有用的。教育的好处显而易见，因为我们可以逐列说明因式分解/高斯消去。供调试使用，这是两个示例。

在Inconsistent results between LU decomposition in R and Python中，询问为什么R和Python中的LU分解会给出不同的结果。我们可以清楚地看到，两个软件都返回相同的第一个枢轴和第二个枢轴，而不是第三个枢轴。因此，当分解进行到第3行/第3列时，一定会有一些有趣的事情。如果我们可以检索该临时结果进行调查，那就太好了。

在Can I stably invert a Vandermonde matrix with many small values in R?中，LU分解对于这种类型的矩阵是不稳定的。在我的回答中，给出了一个3 x 3的矩阵作为示例。我希望solve生成一个错误消息，提示U[3, 3] = 0，但是运行solve几次后，我发现solve有时是成功的。因此，对于数值研究，我想知道当分解进行到第二列/第二行时会发生什么。

由于该函数是用纯R代码编写的，因此对于中等到较大的矩阵，预期它会变慢。但是性能不是问题，因为对于教育和调试，我们只使用一个小的矩阵。

dgetrf的简短介绍

LAPACK的dgetrf通过行数据透视表A = PLU计算LU分解。在因式分解退出时，

L是一个单位下三角矩阵，存储在A的下三角部分中；

U是一个上三角矩阵，存储在A的上三角部分中；

P是存储为单独的排列索引向量的行排列矩阵。

除非支点的恰好为零(不达到一定的容差)，否则应进行分解。

我从开始

编写既不具有行透视功能又不具有“暂停/继续”选项的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

与lu包中的Matrix的比较:

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

请注意，返回的排列索引实际上是针对PA = LU(可能是教科书中使用最多的):

all.equal( B[rr$P, ], with(rr, L %*% U) )
#[1] TRUE

要获取LAPACK返回的排列索引，即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的属性并返回。稍后，当我们将此临时结果传递回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

它告诉您什么也无法做，并按原样返回输入。
最后，让我们尝试一个奇异的方阵。

## 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 - 编写一个可跟踪的R函数，该函数模仿LAPACK的dgetrf进行LU分解，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/51687808/