[6303] in s-news-athena
g-inv and zeros
daemon@ATHENA.MIT.EDU (Pierre Duchesne)
Sun Jan 29 23:01:38 1995
From: Pierre Duchesne <duchesne@STAT.UMontreal.CA>
To: S-news@utstat.toronto.edu
Date: Sun, 29 Jan 1995 22:49:18 -0500 (EST)
Hi,
First of all, can S perform generalized inverse? In the following, i will
suppose that S can't.
Suppose that you want find the Moore-Penrose inverse for that matrix:
> mat
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
Clearly that matrix doesn't have any 'real' inverse:
> solve(mat)
Error in solve.qr(a): apparently singular matrix
Dumped
So there was my first program (i suppose A square... it's enough for me):
> penrose
function(A)
{
U <- svd(A)$u
V <- svd(A)$v
D.nonnul <- svd(A)$d[svd(A)$d!=0]
taille <- length(D.nonnul)
taille.A <- length(A[1, ])
taille2 <- taille.A - taille
invD <- 1/D.nonnul
D <- diag(c(invD, rep(0, taille2)))
G <- V %*% D %*% t(U)
G
}
However, there's a serious bug in that program. For exemple, with our
matrix 'mat' we find that:
> penrose(mat)
[,1] [,2] [,3]
[1,] 3.006737e+14 -6.013474e+14 3.006737e+14
[2,] -6.013474e+14 1.202695e+15 -6.013474e+14
[3,] 3.006737e+14 -6.013474e+14 3.006737e+14
>
Or we know that mat %*% g.mat %*% mat = mat and we find:
> mat %*% g.mat %*% mat
[,1] [,2] [,3]
[1,] -0.25 2 4.25
[2,] 1.00 4 7.00
[3,] 4.50 12 19.50
If we look at the singular value decomposition of 'mat', we find that:
> svd(mat)
$d:
[1] 1.684810e+01 1.068370e+00 5.543107e-16
^^^^^^^^^^^^
So we have a latent value = 0, but the test doesn't pass well in the program...
My first modification is to include a tolerance value (question:
when can you tell that a numeric value in a program is 0???) of 1e-14
> penrose2
function(A)
{
tolerance <- 1e-14
U <- svd(A)$u
V <- svd(A)$v
D.nonnul <- svd(A)$d[svd(A)$d > tolerance]
taille <- length(D.nonnul)
taille.A <- length(A[1, ])
taille2 <- taille.A - taille
invD <- 1/D.nonnul
D <- diag(c(invD, rep(0, taille2)))
G <- V %*% D %*% t(U)
G
}
So now my g-inv of 'mat' become:
> penrose2(mat)
[,1] [,2] [,3]
[1,] -0.6388889 -5.555556e-02 0.5277778
[2,] -0.1666667 4.857226e-17 0.1666667
[3,] 0.3055556 5.555556e-02 -0.1944444
and it's suppose to be the good answer.
But some simulations show that maybe 'sometimes', that program don't give
the right answer...
So if someone, somewhere, has a better program that mine, can you explain
me what you do and send me the code? I will be VERY very appreciated...
Thanks in advance
Pierre Duchesne
duchesnp@nord.stat.umontreal.ca
