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

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

转载作者：行者123 更新时间：2023-12-04 18:56:23

R 核心中没有 LU 分解函数。虽然这样的分解是solve的一步，它没有明确作为独立功能可用。我们可以为此编写一个 R 函数吗？它需要模仿 LAPACK 例程 dgetrf . Matrix包裹有一个 lu function这很好，但如果我们能写一个 会更好可追踪 R函数，可以

将矩阵分解到某一列/行并返回中间结果；

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

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

在 Inconsistent results between LU decomposition in R and Python ，有人问为什么 R 和 Python 中的 LU 分解会给出不同的结果。我们可以清楚地看到，两个软件都返回相同的第 1 个枢轴和第 2 个枢轴，但不是第 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 使用行旋转计算 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,] . | . .

最佳答案

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

带行旋转
这不是太难。
假设 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 factor 不是上三角。这是一个帮助函数来提取它。

## 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/

文章推荐： openssl - 如何创建自己的扩展验证证书以显示绿色条？

文章推荐： bokeh - 如何在 Bokeh 0.13 中制作和弦图？

r - 通过 QR 分解、SVD(和 Cholesky 分解？)计算投影/帽子矩阵
我正在尝试在 R 中计算任意 N x J 矩阵 S 的投影矩阵 P: P = S (S'S) ^ -1 S' 我一直在尝试使用以下函数来执行此操作: P 概述 solve 基于一般方阵的 LU 分解
php - 具有不规则空格和制表符的文件按列拆分/分解
所以我有一个包含数千行的非常旧的文件(我猜是手工生成的)，我正试图将它们移动到一个 rdb 中，但是这些行没有转换为列的格式/模式。例如，文件中的行如下所示: blah blahsdfas
Django:分解 View
这实际上只是一个“最佳实践”问题...... 我发现在开发应用程序时，我经常会得到很多 View 。将这些 View 分解为几个 View 文件是常见的做法吗？换句话说......而不只是有view
r - 分解 `…`参数并分配给多个函数
使用以下函数foo()作为简单示例，如果可能的话，我想将...中给出的值分配给两个不同的函数。 foo args(mapply) function (FUN, ..., MoreArgs = NUL
Groovy 在列表上的不一致解构/分解？
正面案例:可以进入列表 groovy> println GroovySystem.version groovy> final data1 = [[99,2] , [100,4]] groovy> d
algorithm - 数学:分解
省略素数计算方法和因式分解方法的详细信息。为什么要进行因式分解？它的应用是什么？最佳答案哇，这个线程里有这么多争斗。具有讽刺意味的是，这个问题有一个主要的有效答案。因式分解实际上在加密/解
c++ - 分解/重构程序
术语“分解不良”和“重构”程序是什么意思？你能举一个简单的例子来理解基本的区别吗？最佳答案重构是一种通用技术，可以指代许多任务。它通常意味着清理代码、去除冗余、提高代码质量和可读性。分解不良代码
c++ - 分解/分解函数的函数
我以前有，here ，表明 C++ 函数不容易在汇编中表示。现在我有兴趣以一种或另一种方式阅读它们，因为 Callgrind 是 Valgrind 的一部分，在组装时显示它们已损坏。所以我想要么破坏
python - 分解 with 语句
最初，我一直在打开并同时阅读两个文件，内容如下: with open(file1, 'r') as R1: with open(file2, 'r') as R2: ### m
python - Beautifulsoup 分解()
我正在尝试摆脱标签和标签内的内容使用 beatifulsoup。我去看了文档，似乎是一个非常简单的调用函数。有关该功能的更多信息是 here .这是我到目前为止解析的 html 页面的内容...
c# - 分解 float
给定一个 float ，我想将它分成几个部分的总和，每个部分都有给定的位数。例如，给定 3.1415926535 并要求将其分成以 10 为基数的部分，每部分 4 位数字，它将返回 3.141 + 5
jsf - 分解 EAR 文件
我的 JSF 项目被部署为一个 EAR 文件。它还包括一些 war 文件。我需要 EAR 的分解版本(包括分解的内部 WAR)。有什么工具可以做到吗？最佳答案以编程方式还是手动？ EAR 和 W
r - 带有行枢轴的 LU 分解
以下函数不使用行透视进行 LU 分解。 R 中是否有一个现有的函数可以使用行数据进行 LU 分解？ > require(Matrix) > expand(lu(matrix(rnorm(16),4,4
r - 分解 Shiny 代码的最佳实践
关闭。这个问题是opinion-based .它目前不接受答案。想改进这个问题？更新问题，以便 editing this post 提供事实和引用来回答它. 7年前关闭。 Improve this
r - 如何拆分(分解)R中的日期？
我正在使用登记数据进行病假研究。从登记册上，我只得到了每个人的病假开始日期和结束日期。但日期并没有逐年分割。例如，对于人 A，只有开始日期 (1-may-2016) 和结束日期 (14-feb-201
R 高斯消除和 qr 分解
我发现以下 R 代码使用 qr 因式分解无法恢复原始矩阵。我不明白为什么。 a <- matrix(runif(180),ncol=6) a[,c(2,4)] <- 0 b <- qr(a) d <-
r - 用于异常检测的具有缺失值的时间序列的 STL 分解
我正在尝试检测气候数据时间序列中的异常值，其中一些缺失的观测值。在网上搜索我发现了许多可用的方法。其中，STL 分解似乎很有吸引力，因为它去除了趋势和季节性成分并研究了其余部分。阅读 STL: A S
javascript 分解 VIN 号码
我想使用 javascript 分解数组中的 VIN，可能使用正则表达式，然后使用某种循环... 以下是读取 VIN 的方法: http://forum.cardekho.com/topic/600-
scala - 分解 Spark 数据框中的嵌套结构
我正在研究 Databricks 示例。数据框的架构如下所示: > parquetDF.printSchema root |-- department: struct (nullable = true
javascript - 分解 JS 函数
我正在尝试简化我的代码并将其分解为多个文件。例如，我设法做到了: socket.once("disconnect", disconnectSocket); 然后有一个名为 disconnectSock

行者123

个人简介

我是一名优秀的程序员,十分优秀！

作者热门文章

滴滴打车优惠券免费领取

全站热门文章

首页

博学

6Ren·AI

商城

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