[R] extracting the diagonal of an inverse matrix

[William Dunlap]

Wrap the logical
   seq_len(nrow(Asp)) == i
with as.numeric()
   as.numeric(seq_len(nrow(Asp)) == i)
and it will work (slowly).

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


> -----Original Message-----
> From: Camarda, Carlo Giovanni [mailto:Camarda at demogr.mpg.de]
> Sent: Tuesday, April 23, 2013 2:46 AM
> To: William Dunlap; Greg Snow
> Cc: r-help at stat.math.ethz.ch
> Subject: RE: [R] extracting the diagonal of an inverse matrix
> 
> 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 at tibco.com]
> Sent: Monday, April 22, 2013 6:25 PM
> To: Camarda, Carlo Giovanni; Greg Snow
> Cc: r-help at 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 at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> > Of Camarda, Carlo Giovanni
> > Sent: Monday, April 22, 2013 1:51 AM
> > To: Greg Snow
> > Cc: r-help at 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 at gmail.com]
> > Sent: Saturday, April 20, 2013 12:16 AM
> > To: Camarda, Carlo Giovanni
> > Cc: r-help at 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 at demogr.mpg.de<mailto:Camarda at 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 at 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.