updated on this question: I have closed this question and I will post a new question focus on R package Rmpfr.
To conclude this question and to help others, I will post my codes of the inverse of a Vandermonde Matrix from its explicit inverse formula. The generation terms are the x's in [here]1.
I am not a skilled programmer. Therefore I don't expect my codes to be the most efficient one. I post the codes here because it is better than nothing.
library(gtools)
#input is the generation vector of terms of Vandermonde matrix.
FMinv <- function(base){
n=length(base)
inv=matrix(nrow=n,ncol=n)
for (i in 1:n){
for (j in 1:n){
if(j<n){
a=as.matrix(combinations(n,n-j,repeats.allowed = F))
arow.tmp=nrow(a) #this is in fact a[,1]
b=which(a==i)%%length(a[,1])
nrowdel=length(b)
b=replace(b,b==0,length(a[,1]))
a=a[-b,]
if(arow.tmp-nrowdel>1){
a=as.matrix(a)
nrowa=nrow(a)
prod=vector()
for(k in 1:nrowa){
prod[k]=prod(base[a[k,]])
}
num=sum(prod)
}
if(arow.tmp-nrowdel==1){
num=prod(base[a])
}
den=base[i]*prod(base[-i]-base[i])
inv[i,j]=(-1)^(j-1)*num/den
}
if(j==n){
inv[i,j]=1/(base[i]*prod(base[i]-base[-i]))
}
}
}
return(inv)
}
older version
I found an explicit inversion formula.
For a 3 x 3 matrix, the formula will work. Before I try to look further into my situation, can someone help me on how to program this? Calculating it term by term is tiresome.
a <- (1 / 10) ^ c(0, 20, 30)
V <- t(outer(a, 1:(length(a)), "^"))
Vinv = matrix(nrow=3, ncol=3)
Vinv[1,1] = a[2]*a[3]/(a[1]*(a[2]-a[1])*(a[3]-a[1]))
Vinv[1,2] = -(a[2]+a[3])/(a[1]*(a[2]-a[1])*(a[3]-a[1]))
Vinv[1,3] = 1/(a[1]*(a[1]-a[2])*(a[1]-a[3]))
Vinv[2,1] = a[1]*a[3]/(a[2]*(a[1]-a[2])*(a[3]-a[2]))
Vinv[2,2] = -(a[1]+a[3])/(a[2]*(a[1]-a[2])*(a[3]-a[2]))
Vinv[2,3] = 1/(a[2]*(a[2]-a[1])*(a[2]-a[3]))
Vinv[3,1] = a[1]*a[2]/(a[3]*(a[1]-a[3])*(a[2]-a[3]))
Vinv[3,2] = -(a[1]+a[2])/(a[3]*(a[1]-a[3])*(a[2]-a[3]))
Vinv[3,3] = 1/(a[3]*(a[3]-a[1])*(a[3]-a[2]))
det(V %*% Vinv)
#[1] 1
Vi = solve(V)
#Error in solve.default(V) :
# Lapack routine dgesv: system is exactly singular: U[3,3] = 0
old question:
I am using R on some calculations that involves some very small values (like 1e-200).
Suppose I have a full rank n x n
matrix A
. I know I can find the inverse matrix by B <- solve(A, tol = 1e-320)
. Theoretically, det(A %*% B)
should be 1. However, R will give me 0.
I am assuming this is caused by the low precision in %*%
. Is there anyway to reset its precision?
c = 1/(1:50)
A = matrix(nrow = 50, ncol = 50)
for (i in 1:50){
for (j in 1:50){
A[i,j] = c[j] ^ (i-1)
}
}
B = solve(A, tol = 1e-324)
det(A %*% B)
#[1] 0
update on the explicit inverse formula of a Vandermonde matrix
The explicit inverse formula is attractive from a pure math point of view, but is not more stable than LU factorization. You still need to perform floating-point addition, so you still lose significant digits. With the 3 x 3 example matrix in my original answer, the formula needs to compute on the denominator
And for
e1 = 1e-20
ande2 = 1e-30
, numerical computation of1 - e1
and1 - e2
are just 1.You updated your question applying this formula to another 3 x 3 matrix, but have you not realized that numerical computation of, say
are all wrong? So the computed inverse
Vinv
is not analytically correct at all!In fact,
V %*% Vinv
can not be correctly computed, either.Take the
(1, 2)
element for example, you need to sum those values:Looks like
sum(x)
should be 0, butThere is really no point doing all sorts of calculations when they are out of the precision of a machine as the result can just be arbitrary.
original answer on LU factorization of a Vandermonde matrix
You have a Vandermonde matrix, an infamous matrix of ill-conditioning. Let me just quickly generate a small such matrix.
Run this
solve
a few times, you may get two different errors:Now set
tol = 0
,Run this a few times, it might give you the following
Vi
,or it might give you
solve
has no numerically stable behavior here. Even if you are not stopped by an error, the resultingVi
is just numerically wrong, asVi[1, ]
is zero.You get this, because you force a machine to do something impossible for it to do. Inverting the above
V
is not do-able in finite-precision floating-point computations.Appendix: Markdown (MathJax support needed) for the picture: