递归反转线性系统 - 得到巨大的错误（精度）

Question

我有一个 Ax =b 类型的线性系统 - 其中 A 是一个上三角矩阵。 A的结构定义如下：

    comp.Amat <- function(i,j,prob) ifelse(i > j, 0, dbinom(x=i, size=j, prob=prob))

    prob <- 1/4
    A <- outer(1:50, 1:50 , FUN=function(r,c) comp.Amat(r,c,prob) )

A 中的条目是二项式概率 - 问题是当 A 的大小增长时，对角线条目快速接近 0。

如果我们也将向量 b 定义如下：

    b <- seq(1,50,1);

然后求解(a=A,b=b) - 给出错误：

    "  system is computationally singular: reciprocal condition number = 1.07584e-64" 

这是有道理的，因为对角线元素几乎是 0，所以矩阵变得不可逆。

作为一种变通方法，我编写了以下递归函数 - 它开始计算最后一个对角线条目的值，然后替换前几行中的该值。由于矩阵中的每个条目都是 dbinom(j,i, prob) for j=>i ：我可以通过这种方式获得解决方案。

    solve.for.x.custom <- function(A, b, prob)
    {

      n =length(A[1,])
      m =length(A[,1])

      x = seq(1,n, 1);
      x[x> 0] = -1000;

      calc.inv.Aii <- function(i,j, prob)
      {
        res = (1 / (prob*(1-prob)))^i;
        return(res);


      }

      for (i in m:1 )
      {

        if(i ==m)
        {
  rhs =0;

        }else
        {
          rhs=0;
          for(j in m:(i+1))
          {
            rhs =  dbinom(x=i,size=j,prob=prob)*x[j] + rhs;
          }

        }

        x[i] = (b[i] - rhs)*calc.inv.Aii(i,i);

      }
      print(x)
      return(x)

    }

我的问题是 - 当我将此解决方案 x' 乘以矩阵 A 时，误差 (Ax'- b) 很大。由于我有一个分析解决方案（x_i 中的每个条目都可以描述为二项式概率乘以先前的值） - 我应该得到的错误是每行中的 0-。

由于这些问题，我看到 (1 / (1/a)) 可能不等于 a。但是，当前的错误确实很大（-1.13817489781529e+168）。

    x_prime=solve.for.x.custom(A, b, prob)
    A%*%x_prime - b
    #output
                    [,1]
     [1,] -1.13817489781529e+168
     [2,]  2.11872209742428e+167
     [3,] -1.58403954589004e+166
     [4,]  6.52328959209082e+164
     [5,] -1.69562573261261e+163
     [6,]  3.00614551450976e+161
    ***
    [49,]  -7.58010305220250e+08
    [50,]   9.65162608741321e+03

如果您能推荐任何建议或有效的方法，我将不胜感激。我将 A 和 b 的大小设为 50 - 但我也打算增大它们，因此在这种情况下，错误也会增加。

Answer 1

如果您的矩阵A 是上三角矩阵，您可能希望使用backsolve(A, b) 而不是solve(A, b)。

---

您可以在 R 中使用Rmpfr 进行任意精度，这需要编写backsolve 的兼容版本。通过break下面的代码我们可以得到

> print(max(abs(b - .b)), digits=5)
1 'mpfr' number of precision  1024   bits 
[1] 2.9686e-267

但有一个重要的警告：A 中的值可能不够准确，因为它们来自dbinom，而不是使用mpfr 对象。根据您的最终目标，您可能需要使用Rmpfr 编写自己的dbinom 版本。

---

library(Rmpfr)

logcomp.Amat <- function(i,j,prob) ifelse(i > j, -Inf, dbinom(x=i, size=j, prob=prob, log=TRUE))

nbits <- 1024

.backsolve <- function(A, b) {

    n <- length(b)
    x <- mpfr(numeric(n), nbits)

    for(i in rev(seq_len(n))) {

        known <- i + seq_len(n - i)
        z <- if(length(known) > 0) sum(A[i,known] * x[known]) else 0

        x[i] <- (b[i] - z) / A[i,i]
    }

    return(x)
}

logA <- outer(1:50, 1:50, logcomp.Amat, prob=1/4)
b <- 1:50

A <- exp(mpfr(logA, nbits))
b <- mpfr(b, nbits)

x <- .backsolve(A, b)

.b <- as.vector(A %*% x)