[R] Optimization algorithm to be applied to S4 classes - specifically sparse matrices

Ravi Varadhan RVaradhan at jhmi.edu
Fri May 15 20:48:13 CEST 2009


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: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html



----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of spencerg
Sent: Friday, May 15, 2009 2:22 PM
To: Avraham.Adler at guycarp.com
Cc: r-help at 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.Adler at 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

>              <bates at stat.wisc.

>              edu>                                                       To

>              Sent by:                  Avraham.Adler at guycarp.com

>              dmbates at gmail.com                                          cc

>                                        r-help at r-project.org

>                                                                    Subject

>              05/15/2009 11:57          Re: [R] Optimization algorithm to

>              AM                        be applied to S4 classes -

>                                        specifically sparse matrices

>

>

>

>

>

>

>
>
>
>
> On Wed, May 13, 2009 at 5:21 PM,  <Avraham.Adler at 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 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.
>>
>>     
>
>
>
>

______________________________________________
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.



More information about the R-help mailing list