Re: [R] Matrix inversion-different answers from LAPACK and LINPACK
that this is a good way to calculate the result - in fact it is usually a very bad way. whenever a user asks for solve(x) I was corresponding with Tim Davis, an renowned numerical analyst who wrote the sparse matrix software used in the Matrix package, and he was horrified that we even allow the one-argument form of the solve function. He said, But people will do very stupid things if you provide them with an easy way of asking for a matrix inverse and I said, Yep. I would amend fortune(rethink) If the answer is parse() you should usually rethink the question. -- Thomas Lumley R-help (February 2005) to say parse() or solve(x) albyn On Wed, Jun 17, 2009 at 11:37:48AM -0400, avraham.ad...@guycarp.com wrote: Hello. I am trying to invert a matrix, and I am finding that I can get different answers depending on whether I set LAPACK true or false using qr. I had understood that LAPACK is, in general more robust and faster than LINPACK, so I am confused as to why I am getting what seems to be invalid answers. The matrix is ostensibly the Hessian for a function I am optimizing. I want to get the parameter correlations, so I need to invert the matrix. There are times where the standard solve(X) returns an error, but solve(qr(X, LAPACK=TRUE)) returns values. However, there are times, where the latter returns what seems to be bizarre results. For example, an example matrix is PLLH (Pareto LogLikelihood Hessian) alpha theta alpha 1144.6262175141619082 -0.01290205604828788 theta -0.012902056048 0.00155437688485563 Running plain solve (PLLH) or solve (qr(PLLH)) returns: [,1] [,2] alpha 0.0137466171688024 1141.53956787721 theta 1141.5395678772131305 101228592.41439932585 However, running solve(qr(PLLH, LAPACK=TRUE)) returns: [,1] [,2] [1,] 0.0137466171688024 0.0137466171688024 [2,] 1141.5395678772131305 1141.5395678772131305 which seems wrong as the original matrix had identical entries on the off diagonals. I am neither a programmer nor an expert in matrix calculus, so I do not understand why I should be getting different answers using different libraries to perform the ostensibly same function. I understand the extremely small value of d²X/d(theta)² (PLLH[2,2]) may be contributing to the error, but I would appreciate confirmation, or correction, from the experts on this list. Thank you very much, --Avraham Adler PS: For completeness, the QR decompositions per R under LINPACK and LAPACK are shown below: qr(PLLH) $qr alpha theta alpha -1144.6262175869414932095 0.0129020653695122277 theta -0.112768491646264 0.987863187747112 $rank [1] 2 $qraux [1] 1.993641619511209 0.987863187747112 $pivot [1] 1 2 attr(,class) [1] qr qr(PLLH, LAPACK=TRUE) $qr alpha theta alpha -1144.62621758694149320945 0.0129020653695122277 theta -0.0563842458249248 0.987863187747112 $rank [1] 2 $qraux [1] 1.993642 0.00 $pivot [1] 1 2 attr(,useLAPACK) [1] TRUE attr(,class) [1] qr __ R-help@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. __ R-help@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. __ R-help@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. __ R-help@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.
[R] Matrix inversion-different answers from LAPACK and LINPACK
Hello. I am trying to invert a matrix, and I am finding that I can get different answers depending on whether I set LAPACK true or false using qr. I had understood that LAPACK is, in general more robust and faster than LINPACK, so I am confused as to why I am getting what seems to be invalid answers. The matrix is ostensibly the Hessian for a function I am optimizing. I want to get the parameter correlations, so I need to invert the matrix. There are times where the standard solve(X) returns an error, but solve(qr(X, LAPACK=TRUE)) returns values. However, there are times, where the latter returns what seems to be bizarre results. For example, an example matrix is PLLH (Pareto LogLikelihood Hessian) alphatheta alpha 1144.6262175141619082 -0.01290205604828788 theta -0.012902056048 0.00155437688485563 Running plain solve (PLLH) or solve (qr(PLLH)) returns: [,1] [,2] alpha0.0137466171688024 1141.53956787721 theta 1141.5395678772131305 101228592.41439932585 However, running solve(qr(PLLH, LAPACK=TRUE)) returns: [,1] [,2] [1,]0.01374661716880240.0137466171688024 [2,] 1141.5395678772131305 1141.5395678772131305 which seems wrong as the original matrix had identical entries on the off diagonals. I am neither a programmer nor an expert in matrix calculus, so I do not understand why I should be getting different answers using different libraries to perform the ostensibly same function. I understand the extremely small value of d²X/d(theta)² (PLLH[2,2]) may be contributing to the error, but I would appreciate confirmation, or correction, from the experts on this list. Thank you very much, --Avraham Adler PS: For completeness, the QR decompositions per R under LINPACK and LAPACK are shown below: qr(PLLH) $qr alpha theta alpha -1144.6262175869414932095 0.0129020653695122277 theta-0.112768491646264 0.987863187747112 $rank [1] 2 $qraux [1] 1.993641619511209 0.987863187747112 $pivot [1] 1 2 attr(,class) [1] qr qr(PLLH, LAPACK=TRUE) $qr alpha theta alpha -1144.62621758694149320945 0.0129020653695122277 theta-0.0563842458249248 0.987863187747112 $rank [1] 2 $qraux [1] 1.993642 0.00 $pivot [1] 1 2 attr(,useLAPACK) [1] TRUE attr(,class) [1] qr __ R-help@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.
Re: [R] Matrix inversion-different answers from LAPACK and LINPACK
, but I would appreciate confirmation, or correction, from the experts on this list. Thank you very much, --Avraham Adler PS: For completeness, the QR decompositions per R under LINPACK and LAPACK are shown below: qr(PLLH) $qr alpha theta alpha -1144.6262175869414932095 0.0129020653695122277 theta-0.112768491646264 0.987863187747112 $rank [1] 2 $qraux [1] 1.993641619511209 0.987863187747112 $pivot [1] 1 2 attr(,class) [1] qr qr(PLLH, LAPACK=TRUE) $qr alpha theta alpha -1144.62621758694149320945 0.0129020653695122277 theta-0.0563842458249248 0.987863187747112 $rank [1] 2 $qraux [1] 1.993642 0.00 $pivot [1] 1 2 attr(,useLAPACK) [1] TRUE attr(,class) [1] qr __ R-help@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. __ R-help@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.
Re: [R] Matrix inversion-different answers from LAPACK and LINPACK
-boun...@r-project.org] On Behalf Of avraham.ad...@guycarp.com Sent: Wednesday, June 17, 2009 11:38 AM To: r-help@r-project.org Subject: [R] Matrix inversion-different answers from LAPACK and LINPACK Hello. I am trying to invert a matrix, and I am finding that I can get different answers depending on whether I set LAPACK true or false using qr. I had understood that LAPACK is, in general more robust and faster than LINPACK, so I am confused as to why I am getting what seems to be invalid answers. The matrix is ostensibly the Hessian for a function I am optimizing. I want to get the parameter correlations, so I need to invert the matrix. There are times where the standard solve(X) returns an error, but solve(qr(X, LAPACK=TRUE)) returns values. However, there are times, where the latter returns what seems to be bizarre results. For example, an example matrix is PLLH (Pareto LogLikelihood Hessian) alphatheta alpha 1144.6262175141619082 -0.01290205604828788 theta -0.012902056048 0.00155437688485563 Running plain solve (PLLH) or solve (qr(PLLH)) returns: [,1] [,2] alpha0.0137466171688024 1141.53956787721 theta 1141.5395678772131305 101228592.41439932585 However, running solve(qr(PLLH, LAPACK=TRUE)) returns: [,1] [,2] [1,]0.01374661716880240.0137466171688024 [2,] 1141.5395678772131305 1141.5395678772131305 which seems wrong as the original matrix had identical entries on the off diagonals. I am neither a programmer nor an expert in matrix calculus, so I do not understand why I should be getting different answers using different libraries to perform the ostensibly same function. I understand the extremely small value of d²X/d(theta)² (PLLH[2,2]) may be contributing to the error, but I would appreciate confirmation, or correction, from the experts on this list. Thank you very much, --Avraham Adler PS: For completeness, the QR decompositions per R under LINPACK and LAPACK are shown below: qr(PLLH) $qr alpha theta alpha -1144.6262175869414932095 0.0129020653695122277 theta-0.112768491646264 0.987863187747112 $rank [1] 2 $qraux [1] 1.993641619511209 0.987863187747112 $pivot [1] 1 2 attr(,class) [1] qr qr(PLLH, LAPACK=TRUE) $qr alpha theta alpha -1144.62621758694149320945 0.0129020653695122277 theta-0.0563842458249248 0.987863187747112 $rank [1] 2 $qraux [1] 1.993642 0.00 $pivot [1] 1 2 attr(,useLAPACK) [1] TRUE attr(,class) [1] qr __ R-help@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. __ R-help@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.
Re: [R] Matrix inversion-different answers from LAPACK and LINPACK
diagonals. I am neither a programmer nor an expert in matrix calculus, so I do not understand why I should be getting different answers using different libraries to perform the ostensibly same function. I understand the extremely small value of d²X/d(theta)² (PLLH[2,2]) may be contributing to the error, but I would appreciate confirmation, or correction, from the experts on this list. Thank you very much, --Avraham Adler PS: For completeness, the QR decompositions per R under LINPACK and LAPACK are shown below: qr(PLLH) $qr alpha theta alpha -1144.6262175869414932095 0.0129020653695122277 theta -0.112768491646264 0.987863187747112 $rank [1] 2 $qraux [1] 1.993641619511209 0.987863187747112 $pivot [1] 1 2 attr(,class) [1] qr qr(PLLH, LAPACK=TRUE) $qr alpha theta alpha -1144.62621758694149320945 0.0129020653695122277 theta -0.0563842458249248 0.987863187747112 $rank [1] 2 $qraux [1] 1.993642 0.00 $pivot [1] 1 2 attr(,useLAPACK) [1] TRUE attr(,class) [1] qr __ R-help@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. __ R-help@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. __ R-help@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.
Re: [R] Optimization algorithm to be applied to S4 classes - specifically sparse matrices
Thank you both very much for your replies. What makes this a little less straightforward, at least to me, is that there needs to be constraints on the solved parameters. They most certainly need to be positive and there may be an upper limit as well. The true best linear fit would have negative entries for some of the parameters. Originally, I was using the L-BFGS-B method of optim which both allows for box constraints and has the limited memory advantage useful when dealing with large matrices. Having the analytic gradient, I thought of using BFGS and having a statement in the function returning Inf for any parameters outside the allowable constraints. I do /not/ know how to apply parameter constraints when using linear models. I looked around at the various manuals and help features, and outside of package glmc I did not find anything I could use. Perhaps I overlooked something. If there is something I missed, please let me know. If there truly is no standard optimization routine that works on sparse matrices, my next step may be to use the normal equations to shrink the size of the matrix, recast it as a dense matrix (it would only be 1173x1173 then) and then hand it off to optim. Any further suggestions or corrections would be very much appreciated. Thank you, --Avraham Adler Douglas Bates ba...@stat.wisc. edu To Sent by: avraham.ad...@guycarp.com dmba...@gmail.com cc r-help@r-project.org Subject 05/15/2009 11:57 Re: [R] Optimization algorithm to AMbe applied to S4 classes - specifically sparse matrices On Wed, May 13, 2009 at 5:21 PM, avraham.ad...@guycarp.com wrote: Hello. I am trying to optimize a set of parameters using /optim/ in which the actual function to be minimized contains matrix multiplication and is of the form: SUM ((A%*%X - B)^2) where A is a matrix and X and B are vectors, with X as parameter vector. As Spencer Graves pointed out, what you are describing here is a linear least squares problem, which has a direct (i.e. non-iterative) solution. A comparison of the speed of various ways of solving such a system is given in one of the vignettes in the Matrix package. This has worked well so far. Recently, I was given a data set A of size 360440 x 1173, which could not be handled as a normal matrix. I brought it into 'R' as a sparse matrix (dgCMatrix - using sparseMatrix from the Matrix package), and the formulæ and gradient work, but /optim/ returns an error of the form no method for coercing this S4 class to a vector. If you just want the least squares solution X then X - solve(crossprod(A), crossprod(A, B)) will likely be the fastest method where A is the sparse matrix. I do feel obligated to point out that the least squares solution for such large systems is rarely a sensible solution to the underlying problem. If you have over 1000 columns in A and it is very sparse then likely at least parts of A are based on indicator columns for a categorical variable. In such situations a model with random effects for the category is often preferable to the fixed-effects model you are fitting. After briefly looking into methods and classes, I realize I am in way over my head. Is there any way I could use /optim/ or another optimization algorithm, on sparse matrices? Thank you very much, --Avraham Adler __ R-help@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. __ R-help@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.
[R] Optimization algorithm to be applied to S4 classes - specifically sparse matrices
Hello. I am trying to optimize a set of parameters using /optim/ in which the actual function to be minimized contains matrix multiplication and is of the form: SUM ((A%*%X - B)^2) where A is a matrix and X and B are vectors, with X as parameter vector. This has worked well so far. Recently, I was given a data set A of size 360440 x 1173, which could not be handled as a normal matrix. I brought it into 'R' as a sparse matrix (dgCMatrix - using sparseMatrix from the Matrix package), and the formulæ and gradient work, but /optim/ returns an error of the form no method for coercing this S4 class to a vector. After briefly looking into methods and classes, I realize I am in way over my head. Is there any way I could use /optim/ or another optimization algorithm, on sparse matrices? Thank you very much, --Avraham Adler __ R-help@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.
