我可以在 R 中稳定地反转具有许多小值的 Vandermonde 矩阵吗？
Can I stably invert a Vandermonde matrix with many small values in R?
更新了这个问题：
我已经关闭了这个问题，我将 post 一个新问题关注 R 包 Rmpfr。
为了结束这个问题并帮助其他人，我将 post 我的范德蒙德矩阵的逆代码来自其显式逆公式。生成项是 [此处]1.
中的 x
我不是一个熟练的程序员。因此，我不希望我的代码是最有效的代码。我post这里的代码因为它总比没有好。
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)
}
旧版本
我找到了 an explicit inversion formula。
对于 3 x 3 矩阵，该公式有效。在我尝试进一步了解我的情况之前，有人可以帮助我如何编程吗？逐项计算很烦人。
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
老问题：
我在一些涉及一些非常小的值（如 1e-200）的计算中使用 R。
假设我有一个满秩 n x n 矩阵 A。我知道我可以通过 B <- solve(A, tol = 1e-320) 找到逆矩阵。理论上，det(A %*% B)应该是1。但是，R会给我0。
我假设这是由于 %*% 中的低精度造成的。有没有办法重置它的精度？
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
更新 Vandermonde 矩阵的显式逆公式
The explicit inverse formula 从纯数学的角度来看很有吸引力，但并不比 LU 分解更稳定。您仍然需要执行浮点加法，因此您仍然会丢失有效数字。使用我原始答案中的 3 x 3 示例矩阵，公式需要计算分母
1 - e1, 1 - e2, e1 - e2
而对于e1 = 1e-20和e2 = 1e-30，1 - e1和1 - e2的数值计算只是1。
您更新了将此公式应用于另一个 3 x 3 矩阵的问题，但您是否没有意识到
的数值计算
a[1] + a[2], a[2] - a[1], a[3] - a[1]
都错了吗？所以计算出的逆 Vinv 在分析上根本不正确！
其实V %*% Vinv也无法正确计算
round(V %*% Vinv)
# [,1] [,2] [,3]
#[1,] 1 -16384 16384
#[2,] 0 1 0
#[3,] 0 0 1
以 (1, 2) 元素为例，您需要对这些值求和：
x <- V[1, ] * Vinv[, 2]
x
#[1] -1e-20 1e+20 -1e+20
看起来 sum(x) 应该是 0，但是
sum(x)
#[1] -16384
sum(x * 1e-20) * 1e20
#[1] -22204.46
scal <- max(abs(x))
sum(x / scal) * scal
#[1] -11102.23
当超出机器的精度时，进行各种计算真的没有意义，因为结果可能是任意的。
范德蒙矩阵 LU 分解的原始答案
你有一个 Vandermonde matrix，一个臭名昭著的病态矩阵。让我快速生成一个这样的小矩阵。
a <- (1 / 10) ^ c(0, 20, 30)
V <- outer(a, 0:(length(a) - 1), "^")
# [,1] [,2] [,3]
#[1,] 1 1e+00 1e+00
#[2,] 1 1e-20 1e-40
#[3,] 1 1e-30 1e-60
solve(V)
运行这个solve几次，你可能会得到两个不同的错误：
Error in solve.default(V) : system is computationally singular: reciprocal condition number = 1.66667e-61
Error in solve.default(V) : Lapack routine dgesv: system is exactly singular: U[3,3] = 0
现在设置tol = 0,
Vi <- solve(V, tol = 0)
运行这几次，可能会给你下面的Vi,
# [,1] [,2] [,3]
#[1,] 0 0e+00 0e+00
#[2,] 1 1e+60 -1e+60
#[3,] 0 -1e+60 1e+60
或者它可能会给你
Error in solve.default(V) : Lapack routine dgesv: system is exactly singular: U[3,3] = 0
solve 在这里没有数值稳定的行为。即使您没有被错误阻止，结果 Vi 也只是数字错误，因为 Vi[1, ] 为零。
你明白了，因为你强迫机器做一些它不可能做的事情。在有限精度浮点计算中，反转上述 V 是不可行的。
附录： 图片的 Markdown（需要 MathJax 支持）：
Consider the LU of a Vandermonde matrix:
\begin{equation}
V = \begin{pmatrix}
1 & 1 & 1\
1 & e_1 & e_1^2\
1 & e_2 & e_2^2
\end{pmatrix} := LU
\end{equation}
or equivalently a Gaussian elimination $L^{-1}V = U$.
**step 0:** initialize $U = V$, $L^{-1} = I$ and $L$ as an unit lower triangular matrix.
\begin{equation}
U = \begin{pmatrix}
1 & 1 & 1\
1 & e_1 & e_1^2\
1 & e_2 & e_2^2
\end{pmatrix}, \quad L^{-1} = \begin{pmatrix}
1 & 0 & 0\
0 & 1 & 0\
0 & 0 & 1
\end{pmatrix}, \quad L = \begin{pmatrix}
1 & 0 & 0\
\star & 1 & 0\
\star & \star & 1
\end{pmatrix}
\end{equation}
**step 1:** For $U$ and $L^{-1}$, add row 1 $\times\ (-1)$ to row 2 and row 3. For $L$, fill in $L(2,1) = L(3,1) = 1$.
\begin{equation}
U = \begin{pmatrix}
1 & 1 & 1\
0 & e_1 - 1 & e_1^2 - 1\
0 & e_2 - 1 & e_2^2 - 1
\end{pmatrix},\quad L^{-1} = \begin{pmatrix}
1 & 0 & 0\
-1 & 1 & 0\
-1 & 0 & 1
\end{pmatrix},\quad L =\begin{pmatrix}
1 & 0 & 0\
1 & 1 & 0\
1 & \star & 1
\end{pmatrix}
\end{equation}
In finite-precision floating-point computations, for $e_1 < 2.22\times10^{16}$, there is just $e_1 - 1 = -1$. So when $e_1 = 10^{-20}$ and $e_2 = 10 ^{-30}$, these matrices are known to a machine as
\begin{equation}
U = \begin{pmatrix}
1 & 1 & 1\
0 & - 1 & - 1\
0 & - 1 & - 1
\end{pmatrix},\quad L^{-1} = \begin{pmatrix}
1 & 0 & 0\
-1 & 1 & 0\
-1 & 0 & 1
\end{pmatrix},\quad L =\begin{pmatrix}
1 & 0 & 0\
1 & 1 & 0\
1 & \star & 1
\end{pmatrix}
\end{equation}
**step 2:** For $U$ and $L^{-1}$, add row 2 $\times\ (-1)$ to row 3. For $L$, fill in $L(3,2) = 1$.
\begin{equation}
U = \begin{pmatrix}
1 & 1 & 1\
0 & -1 & -1\
0 & 0 & 0
\end{pmatrix}, \quad
L^{-1} = \begin{pmatrix}
1 & 0 & 0\
-1 & 1 & 0\
0 & -1 & 1
\end{pmatrix},\quad L =\begin{pmatrix}
1 & 0 & 0\
1 & 1 & 0\
1 & 1 & 1
\end{pmatrix}
\end{equation}
You will see that
- $U$ is not invertible, so $V^{-1} = U^{-1}L^{-1}$ does not exist.
- $LU = \begin{pmatrix}1 & 1 & 1\1 & 0 & 0\1 & 0 & 0\end{pmatrix}\neq V$
更新了这个问题： 我已经关闭了这个问题，我将 post 一个新问题关注 R 包 Rmpfr。
为了结束这个问题并帮助其他人，我将 post 我的范德蒙德矩阵的逆代码来自其显式逆公式。生成项是 [此处]1.
中的 x我不是一个熟练的程序员。因此，我不希望我的代码是最有效的代码。我post这里的代码因为它总比没有好。
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)
}
旧版本
我找到了 an explicit inversion formula。
对于 3 x 3 矩阵，该公式有效。在我尝试进一步了解我的情况之前，有人可以帮助我如何编程吗？逐项计算很烦人。
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
老问题：
我在一些涉及一些非常小的值（如 1e-200）的计算中使用 R。
假设我有一个满秩 n x n 矩阵 A。我知道我可以通过 B <- solve(A, tol = 1e-320) 找到逆矩阵。理论上，det(A %*% B)应该是1。但是，R会给我0。
我假设这是由于 %*% 中的低精度造成的。有没有办法重置它的精度？
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
更新 Vandermonde 矩阵的显式逆公式
The explicit inverse formula 从纯数学的角度来看很有吸引力，但并不比 LU 分解更稳定。您仍然需要执行浮点加法，因此您仍然会丢失有效数字。使用我原始答案中的 3 x 3 示例矩阵，公式需要计算分母
1 - e1, 1 - e2, e1 - e2
而对于e1 = 1e-20和e2 = 1e-30，1 - e1和1 - e2的数值计算只是1。
您更新了将此公式应用于另一个 3 x 3 矩阵的问题，但您是否没有意识到
的数值计算a[1] + a[2], a[2] - a[1], a[3] - a[1]
都错了吗？所以计算出的逆 Vinv 在分析上根本不正确！
其实V %*% Vinv也无法正确计算
round(V %*% Vinv)
# [,1] [,2] [,3]
#[1,] 1 -16384 16384
#[2,] 0 1 0
#[3,] 0 0 1
以 (1, 2) 元素为例，您需要对这些值求和：
x <- V[1, ] * Vinv[, 2]
x
#[1] -1e-20 1e+20 -1e+20
看起来 sum(x) 应该是 0，但是
sum(x)
#[1] -16384
sum(x * 1e-20) * 1e20
#[1] -22204.46
scal <- max(abs(x))
sum(x / scal) * scal
#[1] -11102.23
当超出机器的精度时，进行各种计算真的没有意义，因为结果可能是任意的。
范德蒙矩阵 LU 分解的原始答案
你有一个 Vandermonde matrix，一个臭名昭著的病态矩阵。让我快速生成一个这样的小矩阵。
a <- (1 / 10) ^ c(0, 20, 30)
V <- outer(a, 0:(length(a) - 1), "^")
# [,1] [,2] [,3]
#[1,] 1 1e+00 1e+00
#[2,] 1 1e-20 1e-40
#[3,] 1 1e-30 1e-60
solve(V)
运行这个solve几次，你可能会得到两个不同的错误：
Error in solve.default(V) : system is computationally singular: reciprocal condition number = 1.66667e-61
Error in solve.default(V) : Lapack routine dgesv: system is exactly singular: U[3,3] = 0
现在设置tol = 0,
Vi <- solve(V, tol = 0)
运行这几次，可能会给你下面的Vi,
# [,1] [,2] [,3]
#[1,] 0 0e+00 0e+00
#[2,] 1 1e+60 -1e+60
#[3,] 0 -1e+60 1e+60
或者它可能会给你
Error in solve.default(V) : Lapack routine dgesv: system is exactly singular: U[3,3] = 0
solve 在这里没有数值稳定的行为。即使您没有被错误阻止，结果 Vi 也只是数字错误，因为 Vi[1, ] 为零。
你明白了，因为你强迫机器做一些它不可能做的事情。在有限精度浮点计算中，反转上述 V 是不可行的。
附录： 图片的 Markdown（需要 MathJax 支持）：
Consider the LU of a Vandermonde matrix:
\begin{equation}
V = \begin{pmatrix}
1 & 1 & 1\
1 & e_1 & e_1^2\
1 & e_2 & e_2^2
\end{pmatrix} := LU
\end{equation}
or equivalently a Gaussian elimination $L^{-1}V = U$.
**step 0:** initialize $U = V$, $L^{-1} = I$ and $L$ as an unit lower triangular matrix.
\begin{equation}
U = \begin{pmatrix}
1 & 1 & 1\
1 & e_1 & e_1^2\
1 & e_2 & e_2^2
\end{pmatrix}, \quad L^{-1} = \begin{pmatrix}
1 & 0 & 0\
0 & 1 & 0\
0 & 0 & 1
\end{pmatrix}, \quad L = \begin{pmatrix}
1 & 0 & 0\
\star & 1 & 0\
\star & \star & 1
\end{pmatrix}
\end{equation}
**step 1:** For $U$ and $L^{-1}$, add row 1 $\times\ (-1)$ to row 2 and row 3. For $L$, fill in $L(2,1) = L(3,1) = 1$.
\begin{equation}
U = \begin{pmatrix}
1 & 1 & 1\
0 & e_1 - 1 & e_1^2 - 1\
0 & e_2 - 1 & e_2^2 - 1
\end{pmatrix},\quad L^{-1} = \begin{pmatrix}
1 & 0 & 0\
-1 & 1 & 0\
-1 & 0 & 1
\end{pmatrix},\quad L =\begin{pmatrix}
1 & 0 & 0\
1 & 1 & 0\
1 & \star & 1
\end{pmatrix}
\end{equation}
In finite-precision floating-point computations, for $e_1 < 2.22\times10^{16}$, there is just $e_1 - 1 = -1$. So when $e_1 = 10^{-20}$ and $e_2 = 10 ^{-30}$, these matrices are known to a machine as
\begin{equation}
U = \begin{pmatrix}
1 & 1 & 1\
0 & - 1 & - 1\
0 & - 1 & - 1
\end{pmatrix},\quad L^{-1} = \begin{pmatrix}
1 & 0 & 0\
-1 & 1 & 0\
-1 & 0 & 1
\end{pmatrix},\quad L =\begin{pmatrix}
1 & 0 & 0\
1 & 1 & 0\
1 & \star & 1
\end{pmatrix}
\end{equation}
**step 2:** For $U$ and $L^{-1}$, add row 2 $\times\ (-1)$ to row 3. For $L$, fill in $L(3,2) = 1$.
\begin{equation}
U = \begin{pmatrix}
1 & 1 & 1\
0 & -1 & -1\
0 & 0 & 0
\end{pmatrix}, \quad
L^{-1} = \begin{pmatrix}
1 & 0 & 0\
-1 & 1 & 0\
0 & -1 & 1
\end{pmatrix},\quad L =\begin{pmatrix}
1 & 0 & 0\
1 & 1 & 0\
1 & 1 & 1
\end{pmatrix}
\end{equation}
You will see that
- $U$ is not invertible, so $V^{-1} = U^{-1}L^{-1}$ does not exist.
- $LU = \begin{pmatrix}1 & 1 & 1\1 & 0 & 0\1 & 0 & 0\end{pmatrix}\neq V$