Re: [R] Optimization algorithm to be 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.
Re: [R] Optimization algorithm to be applied to S4 classes - specifically sparse matrices
Dear Avraham:
For problems with many parameters to estimate, I highly recommend
Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus
(Springer). This book includes numerous examples showing how to use the
nlme package. The value of this book is greatly enhanced by the
availability of script files named, ch01.R, ch02.R, ... ch08.R
showing how to work virtually all the examples in the book. These
script files are available in your local installation of R. To find
them, enter the following at a commands prompt in R:
system.file('scripts', package='nlme')
Hope this helps.
Spencer Graves
##
Dear Doug, et al.:
What would you recommend for analyzing a longitudinal abundance
survey of 22 species, when the species were not selected at random? A
prominent scientist tried to tell me that mixed-effects modeling is
inappropriate in that case because the species were selected
purposefully not at random.
My response is that even in that case, one should still use
mixed-effects modeling, because it will tend to produce more appropriate
estimates for the deviations of individual species from the average of
all species -- potentially much lower variance with slight bias -- than
naive ordinary least squares. The estimated variance components will
not represent the between-species variance for the actual population of
all hypothetical species of the particular type, but will represent the
between-species variability in a hypothetical population from which the
selected species might be considered a random sample.
Best Wishes,
Spencer Graves
p.s. I appreciate very much Doug's comment on this. I thought about
adding something like that to my reply but didn't feel I could afford
the time then.
Douglas Bates wrote:
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-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
Dear Doug, et al.: What would you recommend for analyzing a longitudinal abundance survey of 22 species, when the species were not selected at random? A prominent scientist tried to tell me that mixed-effects modeling is inappropriate in that case because the species were selected purposefully not at random. My response is that even in that case, one should still use mixed-effects modeling, because it will tend to produce more appropriate estimates for the deviations of individual species from the average of all species -- potentially much lower variance with slight bias -- than naive ordinary least squares. The estimated variance components will not represent the between-species variance for the actual population of all hypothetical species of the particular type, but will represent the between-species variability in a hypothetical population from which the selected species might be considered a random sample. Best Wishes, Spencer Graves p.s. I appreciate very much Doug's comment on this. I thought about adding something like that to my reply but didn't feel I could afford the time then. Douglas Bates wrote: 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-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.
Re: [R] Optimization algorithm to be applied to S4 classes - specifically sparse matrices
I suggest you try to translate your constraints into an unconstrained constrained problem using logarithms, then do nonlinear mixed effects modeling as described in chapters 6-8 of Pinheiro and Bates (2000). To do this, I would first start with the simpler linear estimation problem to get starting values for the nonlinear estimation. You should be able to do this using the nlme function in the nlme package. If you have trouble with this, you might consider the nlmer function in the lme4 package. The latter is newer and better in many ways but not as well documented. Hope this helps. Spencer Graves avraham.ad...@guycarp.com wrote: 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/
Re: [R] Optimization algorithm to be applied to S4 classes - specifically sparse matrices
Hi, I think quadratic programming is the way to go. Look at solve.QP or limSolve package. Here is a toy example that I had worked out some time back for a linear least squares problem with simple box constraints: # Problem: minimize ||Ax - y||, subject to low = x = upp require(limSolve) nc - 7 # 7 unknown parameters nr - 20 # 20 equations # Bounds on the parameters: 0 x 1, for all x # set.seed(123) A - matrix(rnorm(nr*nc), nr, nc) x - c(runif(nc-1), 1.5) # Note: the last component is out of bounds! y - A %*% x + rnorm(nr, sd=0.1) qr.solve(A, y) # unconstrained least-squares low - rep(0, nc) # lower bounds upp - rep(1, nc) # upper bounds # Implementing the bounds (there is probably a simpler way to do this) # c1 - matrix(0, nc, nc) diag(c1) - 1 c2 - matrix(0, nc, nc) diag(c2) - -1 cmat - rbind(c1, c2) vec - rep(0, 10) vec[seq(1, 2*nc, by=2)] - 1:nc vec[seq(2, 2*nc, by=2)] - (nc+1):(2*nc) Cmat - rbind(c1, c2)[vec, ] # Constraint matrix G b0 - c(low, -upp)[vec] ans - lsei(A = A, B = y, G = Cmat, H = b0) ans Hope this helps, Ravi. --- Ravi Varadhan, Ph.D. Assistant Professor, The Center on Aging and Health Division of Geriatric Medicine and Gerontology Johns Hopkins University Ph: (410) 502-2619 Fax: (410) 614-9625 Email: rvarad...@jhmi.edu Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html -Original Message- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of spencerg Sent: Friday, May 15, 2009 2:22 PM To: avraham.ad...@guycarp.com Cc: r-help@r-project.org; Douglas Bates Subject: Re: [R] Optimization algorithm to be applied to S4 classes - specifically sparse matrices I suggest you try to translate your constraints into an unconstrained constrained problem using logarithms, then do nonlinear mixed effects modeling as described in chapters 6-8 of Pinheiro and Bates (2000). To do this, I would first start with the simpler linear estimation problem to get starting values for the nonlinear estimation. You should be able to do this using the nlme function in the nlme package. If you have trouble with this, you might consider the nlmer function in the lme4 package. The latter is newer and better in many ways but not as well documented. Hope this helps. Spencer Graves avraham.ad...@guycarp.com wrote: 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
Re: [R] Optimization algorithm to be applied to S4 classes - specifically sparse matrices
Have you considered the following: solve(qr(A), B) I have not tried this with a small toy example, and the qr documentation in the Matrix package seems to suggest it. This solves the optimization problem you mentioned, as noted in http://en.wikipedia.org/wiki/Linear_least_squares;. Another alternative is the biglm function in the package of the same name. I have not tried this either, but it looks like it should work. Hope this helps. Spencer Graves 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. 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. __ 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.
