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