相关疑难解决方法(0)

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] …

我正在尝试使用 R 运行 LU 分解。这是可重现的代码。我不明白为什么我的排列矩阵与解决方案不同。L 和 U 矩阵是正确的。但对于置换矩阵，第 1 行和第 2 行以及第 3 行和第 4 行互换。因此，我没有得到线性方程组的正确解。将不胜感激您的帮助。

install.packages("Matrix")
library(Matrix)
(A <- matrix(c(4, 3, -2, 5, 2, -4, 6, 1, -1, 2, -5, 6, 3, 5, -2, -3), nrow = 4))
(B <- matrix(c(16.9, -14, 25, 9.4), nrow = 4))

luA <- lu(A)
elu <- expand(luA)
(L <- elu$L)
(U <- elu$U)
(P <- elu$P)

(Y <- solve(L) %*% P %*% B)
(X <- solve(U) %*% Y)