0

私は 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);

次に、solve(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 にしましたが、それらも大きくするつもりなので、その場合はエラーも大きくなります。

4

1 に答える 1

2

If your matrix A is upper triangular you probably want to use backsolve(A, b) rather than solve(A, b).

You can do arbitrary precision in R with Rmpfr, which will require writing a compatible version of backsolve. With the code below the break we can get

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

There is one important caveat though: the values in A may not be accurate enough since they come from dbinom rather than using mpfr objeccts. Depending on your end goal, you may need to write your own version of dbinom using Rmpfr.

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)
于 2012-11-10T01:07:21.860 に答える