我有一个 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 - 但我也打算增大它们,因此在这种情况下,错误也会增加。