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

Question

R核心中没有LU分解功能.尽管这种分解是一个步骤solve,但它并未明确地作为独立功能提供.我们可以为此写一个R函数吗？它需要模仿LAPACK例程dgetrf.Matrixpackage有一个很好的lu函数,但如果我们可以编写一个可追踪的 R函数会更好

• 将矩阵分解到某个列/行并返回中间结果;
• 继续从中间结果到另一列/行或到结尾的因子分解.

此功能对于教育和调试目的都很有用.教育的好处是显而易见的,因为我们可以逐列说明分解/高斯消除.对于调试使用,这里有两个例子.

在R和Python中LU分解之间的结果不一致时,人们会问为什么R和Python中的LU分解会产生不同的结果.我们可以清楚地看到,两个软件都返回相同的第一个枢轴和第二个枢轴,但不是第三个.因此,当分解进行到第3行/列时,必定会有一些有趣的东西.如果我们能够检索调查的临时结果,那将是件好事.

在我可以稳定地反转与R中许多小值范德蒙矩阵？对于这种类型的矩阵,LU分解是不稳定的.在我的回答中,给出了一个3 x 3矩阵的例子.我希望solve产生错误抱怨U[3, 3] = 0,但运行solve几次我发现solve有时候会成功.因此,对于数值研究,我想知道当分解进行到第二列/行时会发生什么.

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

---

dgetrf的一点介绍

LAPACK的dgetrf用行旋转计算LU分解:A = PLU.退出分解时,

• 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,] . | . .

Answer 1

让我们逐个添加这些功能.

行旋转

这不是太难.

假设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返回的排列索引,即,in A = PLU,do 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