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

【问题标题】：Write a trackable R function that mimics LAPACK's dgetrf for LU factorization编写一个可跟踪的 R 函数，该函数模仿 LAPACK 的 dgetrf 进行 LU 分解
【发布时间】：2019-01-12 05:50:29
【问题描述】：

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

  • 将矩阵分解到某一列/行并返回中间结果；
  • 继续从中间结果分解到另一列/行或到末尾。

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

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

在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 使用行旋转计算 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,] . | . .

【问题讨论】：

  • pracma 包也有一个 lu 函数，但它不应用旋转。读者可以验证（使用矩阵B）它等同于我的问题中的函数LU。 pracma::lu 在 R 级别使用双循环嵌套编写得很糟糕。 LU 使用单个 R 级循环，因此是更好的实现。

标签： r function matrix matrix-factorization


【解决方案1】：

让我们一一添加这些功能。

行旋转

这并不太难。

假设A 是n x n。初始化一个置换索引向量pivot &lt;- 1:n。在j-th 列，我们扫描A[j:n, j] 以获得最大绝对值。假设它是A[m, j]。如果m &gt; j 我们进行行交换A[m, ] &lt;-&gt; A[j, ]。同时我们做一个排列pivot[j] &lt;-&gt; 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

【讨论】：

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 2020-10-21
    • 1970-01-01
    • 2018-04-29
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多
    热门标签
    Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode