[R] extracting the diagonal of an inverse matrix

[Camarda, Carlo Giovanni]

Thanks: really nice idea. 
Unfortunately the solution you suggest can not be adopted with sparse matrices: solve command with logical is not yet implemented. (Or do I have a too-old R-version?)

## loading package
library(Matrix)
## transform A in sparse matrix
Asp <- Diagonal(9)
pos <- which(A!=0)
Asp[pos] <- A[pos]
## attempt the solution
sapply(seq_len(nrow(Asp)), function(i)solve(Asp, seq_len(nrow(Asp))==i)[i])

## Error: not-yet-implemented method for solve(<dgCMatrix>, <logical>).
##  ->>  Ask the package authors to implement the missing feature.


> version
               _                            
platform       i686-pc-linux-gnu            
arch           i686                         
os             linux-gnu                    
system         i686, linux-gnu              
status                                      
major          2                            
minor          14.1                         
year           2011                         
month          12                           
day            22                           
svn rev        57956                        
language       R                            
version.string R version 2.14.1 (2011-12-22)

________________________________________
From: William Dunlap [wdunlap a tibco.com]
Sent: Monday, April 22, 2013 6:25 PM
To: Camarda, Carlo Giovanni; Greg Snow
Cc: r-help a stat.math.ethz.ch
Subject: RE: [R] extracting the diagonal of an inverse matrix

You could save the space required to store the inverse of A by repeating solve(A,column)
for each column of an identity matrix the size of A, and store just the relevant element
of the result.  E.g.,

> diag(solve(A))
[1] 0.39576620 0.12064108 0.20705171 0.10067892 0.03191425 0.05857497 0.10749426 0.03868911 0.08550176
> sapply(seq_len(nrow(A)), function(i)solve(A, seq_len(nrow(A))==i)[i])
[1] 0.39576620 0.12064108 0.20705171 0.10067892 0.03191425 0.05857497 0.10749426 0.03868911 0.08550176

It may take more time that you would like.

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


> -----Original Message-----
> From: r-help-bounces a r-project.org [mailto:r-help-bounces a r-project.org] On Behalf
> Of Camarda, Carlo Giovanni
> Sent: Monday, April 22, 2013 1:51 AM
> To: Greg Snow
> Cc: r-help a stat.math.ethz.ch
> Subject: Re: [R] extracting the diagonal of an inverse matrix
>
> Dear Greg,
>
> thanks a lot for your prompt reply. I was aware of the link and the associated papers.
> In my previous mail I was wondering whether there was already an R-routine for coping
> with this issue.
>
> I add below a testable example which reproduces my problem.
>
> Thanks again,
> Giancarlo
>
> ## dimensions (my actual dimensions are much larger and I use sparse matrices)
> m <- 3
> n <- 3
> ## building A
> ## diagonal part
> d <- 1:(m*n)
> D <- diag(d)
> ## additional ~dense part (:penalty term)
> P1 <- diff(diag(m), diff=2)
> P2 <- diff(diag(n), diff=2)
> P1a <- kronecker(diag(n), t(P1)%*%P1)
> P2a <- kronecker(t(P2)%*%P2, diag(m))
> P <- 1000 * P1a + 1000 * P2a
> ## final matrix A
> A <- D + P
> ## solving A
> B <- solve(A)
> ## what I just actually need
> diag(B)
>
>
> ________________________________________
> From: Greg Snow [538280 a gmail.com]
> Sent: Saturday, April 20, 2013 12:16 AM
> To: Camarda, Carlo Giovanni
> Cc: r-help a stat.math.ethz.ch
> Subject: Re: [R] extracting the diagonal of an inverse matrix
>
> This link http://math.stackexchange.com/questions/18488/diagonal-of-an-inverse-of-a-
> sparse-matrix might help.  It is about sparse matrices, but the general idea should be able
> to be extended to non-sparse matrices as well.
>
>
> On Fri, Apr 19, 2013 at 8:13 AM, Camarda, Carlo Giovanni
> <Camarda a demogr.mpg.de<mailto:Camarda a demogr.mpg.de>> wrote:
> Dear R-users,
>
> I would like to know whether there is a way to extract a diagonal of an inverse matrix
> without computing the inverse of the matrix itself. The size of my matrices are really
> huge and, also using sparse matrix, computing the inverse leads to storage problems and
> low speed.
>
> In other words, given a square matrix A, I aim to know diag(B), where B=solve(A),
> without computing solve(A).
>
> Accidentally (I do not know whether it helps), I could write the matrix A as follows:
> A <- D + P
> where D is a diagonal matrix.
>
> I read there are methods around, but, before implementing one of them by myself, could
> you please inform whether there is already an R-routine for this issue?
>
> Thanks in advance for the help you could provide,
> Carlo Giovanni Camarda
> ----------
> This mail has been sent through the MPI for Demographic ...{{dropped:19}}
>
> ______________________________________________
> R-help a r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
----------
This mail has been sent through the MPI for Demographic Research.  Should you receive a mail that is apparently from a MPI user without this text displayed, then the address has most likely been faked. If you are uncertain about the validity of this message, please check the mail header or ask your system administrator for assistance.