SCM

SCM Repository

[matrix] Diff of /pkg/TODO
ViewVC logotype

Diff of /pkg/TODO

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1087, Fri Dec 9 20:55:39 2005 UTC revision 2175, Wed Apr 23 11:23:50 2008 UTC
# Line 1  Line 1 
1  - Migration of lmer from the mer representation to the mer2  - Check for DimNames propagation in coercion and other operations.
    representation and the use of the CHOLMOD code for the sparse  
    matrix decomposition.  Some of the things that need to be done.  
   
    - Matrices in the mer2 representation are classed matrices, in the  
    mer representation they were unclassed.  Any parts inside the C  
    code that would access, for example,  
     REAL(GET_SLOT(x, Matrix_RXXSym))  
    need to be modified to access  
     REAL(GET_SLOT(GET_SLOT(x, Matrix_RXXSym), Matrix_xSym))  
    This is especially important for Omega but I think I have done  
    those changes already.  
   
    - The components named *X* in an mer object refer to an augmented  
    design matrix of p+1 columns.  In the mer2 object there are  
    separate slots for rZy and rXy.  The scalar y'y is the first  
    element of devComp.  
   
    - Presently nc is of length nf+1 and the last element is n, the  
    number of observations.  This value should be moved to devComp and  
    nc made of length nf.  
   
    - The slot L is a list of length 1 that contains an ExternalPointer  
    to a cholmod_factor object.  This contains a permutation which is  
    most easily accessible through cholmod_solve(CHOLMOD_P,...) or  
    cholmod_solve(CHOLMOD_Pt,...).  The ZtX, Zty, RZX and rZy slots actually  
    contain P%*%RZX and P%*%rZy  
   
 ------  
 - Sparse matrix methods can now be based on the CHOLMOD package.  We  
    will need to migrate from the current code to CHOLMOD-based code  
    using #ifdef USE_CHOLMOD.  Some of the things to be done  
   
    - Move documentation from subdirectories of src to inst/doc  
    - Write utilities to create a cholmod_sparse pointer from a  
    dgCMatrix or lgCMatrix (or zgCMatrix) object without copying and  
    allocating.  
    - Start adding simple S4 methods (rcond, %*%, +, cbind, t).  
2    
3  - Report the problem in the Linux ldexp manual page.  The second and  - Report the problem in the Linux ldexp manual page.  The second and
4    third calls in the Synopsis should be to ldexpf and ldexpl.    third calls in the Synopsis should be to ldexpf and ldexpl.
5    
 - [,] indexing: for sparse "works", but not yet for negative indices!  
   
 - consider moving alloc3Darray from ./src/Mutils.c to  
   $(RSRC)/src/base/array.c  
   
 -------  
   
6  - provide methods for "dspMatrix" and "dppMatrix"!  - provide methods for "dspMatrix" and "dppMatrix"!
7    
8  - implement (more) methods for supporting "packed" (symmetric / triangular)  - implement (more) methods for supporting "packed" (symmetric / triangular)
# Line 55  Line 11 
11    
12    (have some dtr* <-> dtp*)    (have some dtr* <-> dtp*)
13    
 - implement diagonal Matrix class  "ddiMatrix" etc  
   using constructor function Diagonal() or Diag().  
   
 - rcond() of a singular dpoMatrix gives a LaPack error instead of just 0:  
   MM <- crossprod(M <- Matrix(c(1:4,9:6), 2,4)) ; rcond(MM)  
   ##> Error in rcond(MM) : Lapack routine dpotrf returned error code 4  
   It's .Call("dpoMatrix_rcond") --> set_rcond() in src/dpoMatrix.c  
   and in src/dppMatrix.c similarly.  
   
   Done(2005-10-03): The error message is more helpful now.  
   
 ---  
   
14  - combine the C functions for multiplication by special forms and  - combine the C functions for multiplication by special forms and
15    solution wrt special forms by using a 'right' argument and a    solution wrt special forms by using a 'right' argument and a
16    'classed' argument.    'classed' argument.
17     [done with dgeMatrix_matrix_mm();  not yet for other classes;     [done with dgeMatrix_matrix_mm();  not yet for other classes;
18      and for _crossprod()]      and for _crossprod()]
19    
20  - add more comprehensive examples / tests for Schur decomposition  -----
   
 - arithmetic for sparse matrices:  
              <sparseMatrix>  o  <same-dim-sparseMatrix>  
   should return a sparse matrix  for at least "+" and "*" , also %%,  
   and "/" and "%/%" at least when the RHS is non-zero a scalar.  
   Challenge: nice implementation (``common non-0''; but Tsparse* is not uniq).  
   
 ---  
   
 - Create a Harwell-Boeing version of the matrix mm and the response  
   vector y in inst/external and remove them from the data directory.  
   Modify any examples that use them and modify the Comparisons vignette.  
21    
22  - "Math2" , "Math", "Arith":  - "Math2" , "Math", "Arith":
23     keep triangular and symmetric Matrices when appropriate:     keep triangular and symmetric Matrices when appropriate:
24     particularly desirable for  "Math2": round(), signif()     particularly desirable for  "Math2": round(), signif()
25    
26      For triangular matrices, more specifically make sure the four rules of
27      "triangular matrix algebra" (Golub+Van Loan 1996, 3.1.8, p.93) are
28      fulfilled; now(2008-03-06) ok for Csparse; not yet for <dtr> %*% <dtr>
29    
30  - "d" <-> "l" coercion for all "[TCR]" sparse matrices is really trivial:  - "d" <-> "l" coercion for all "[TCR]" sparse matrices is really trivial:
31    "d" -> "l" : drops the 'x' slot    "d" -> "l" : drops the 'x' slot
32    "l" -> "d" : construct an 'x' slot of all '1'    "l" -> "d" : construct an 'x' slot of all '1'
# Line 104  Line 39 
39    for all  "dsparse*" -> "lsparse*" and vice versa.    for all  "dsparse*" -> "lsparse*" and vice versa.
40    How can one do this {in a documented way} ?    How can one do this {in a documented way} ?
41    
42  - tcrossprod(x, y) : do provide methods for y != NULL  -- at least dummy ones!  - Think of constructing  setAs(...) calls automatically in order to
43      basically enable all ``sensible'' as(fromMatrix, toMatrix)  calls,
44      possibly using canCoerce(.)
45    
46    - setAs(<Mcl>,  "[dln]Matrix") for <Mcl> in {Matrix or denseMatrix + sparseMatrix}
47    
48    - When we have a packed matrix, it's a waste to go through "full" to "sparse":
49      ==> implement
50            setAs("dspMatrix", "sparseMatrix")
51            setAs("dppMatrix", "sparseMatrix")
52            setAs("dtpMatrix", "sparseMatrix")
53      and the same for "lsp" , "ltp"  and  "nsp" , "ntp" !
54    
55    - tcrossprod(x, y) : do provide methods for y != NULL
56      calling Lapack's DGEMM for "dense"
57      [2005-12-xx: done for dgeMatrix at least]
58    
59    - BUGlet:  Shouldn't lose factorization here:
60      h6 <- Hilbert(6); chol(h6) ; str(h6) # has factor
61      str(H6 <- as(h6, "dspMatrix"))       # has lost factor
62      ## and the same in a similar situation involving  "dpo", "dpp"
63    
64    - Factorizations: LU done; also Schur()  for  *sparse*  Matrices.
65    
66    - is.na() method for all our matrices [ ==> which(*, arr.ind=TRUE) might work ]
67    
68    - use  .Call(Csparse_drop, M, tol) in more places,
69      both with 'tol = 0.' to drop "values that happen to be 0" and for
70      zapsmall() methods for Csparse*
71    
72    - implement .Call(Csparse_scale, ....) interfacing to cholmod_scale()
73      in src/CHOLMOD/Include/cholmod_matrixops.h : for another function
74      specifically for multiplying a cholmod_sparse object by a diagonal matrix.
75      Use it in %*% and [t]crossprod methods.
76    
77    - chol() and determinant() should ``work'': proper result or "good" error
78      message.
79    
80    - make sure *all* group methods have (maybe "bail-out") setMethod for "Matrix".
81      e.g. zapsmall(<pMatrix>) fails "badly"
82    
83    - sum(): implement methods which work for *all* our matrices.
84    
85    - Implement  expand(.) for the Cholesky() results
86      "dCHMsimpl" and  "dCHMsuper"  -- currently have no *decent* way to get at
87      the matrix factors of the corresponding matrix factorization !!
88    
89    - rbind2(<sparse>, <dense>) does not work  (e.g. <dgC>, <dge>)
90    
91    - <sparse> %*% <dense>  {also in crossprod/tcrossprod}  currently always
92      returns <dense>, since --> Csparse_dense_prod --> cholmod_sdmult
93      and that does only return dense.
94      When the sparse matrix is very sparse, i.e. has many rows with only zero
95      entries, it would make much sense to return sparse.
96    
97    - sparse-symmetric + diagonal should stay sparse-symmetric
98      (only stays sparse): Matrix(0, 4, 4) + Diagonal(4, 1:4)
99      --> R/diagMatrix.R ('FIXME')
100      but also R/Ops.R  to ensure  sp-sym. + sp-sym. |-> sp-sym.  etc
101    
102    - Diagonal(n) %*% A ---  too slow!! --> ~/R/MM/Pkg-ex/Matrix/diag-Tamas-ex.R
103    
104    - ! <symmetricMatrix>  loses symmetry, both for dense and sparse matrices.
105      !M  where M is "sparseMatrix", currently always gives dense. This only
106      makes sense when M is ``really sparse''.
107    
108    - msy <- as(matrix(c(2:1,1:2),2), "dsyMatrix"); str(msy)
109    
110      shows that the Cholesky factorization is computed ``too quickly''.
111      Can be a big pain for largish matrices, when it is unneeded.
112    
113    - example(Cholesky, echo=FALSE) ; cm <- chol(mtm); str(cm); str(mtm)
114    
115      shows that chol() does not seems to make use of an already
116      present factorization and rather uses one with more '0' in x slot.
117    
118    - diag(m) <- val    currently automatically works via  m[cbind(i,i)] <- val
119      This (`[<-` method) is now "smart" for diagonalMatrix, but needs also to
120      be for triangularMatrix, and probably also "dense*general*Matrix" since the
121      above currently goes via "matrix" and back instead of using the 'x' slot
122      directly; in particular, the triangular* "class property" is lost!
123    
124    - examples for solve( Cholesky(.), b, system = c("A", "LDLt"....))
125      probably rather in man/CHMfactor-class.Rd than man/Cholesky.Rd
126    
127    - LDL(<CHMsimpl>) looks relatively easy; via  "tCsparse_diag()"
128       {diagonal entries of *triangular* Csparse}
129      --> see comment in determinant(<dsC>) in R/dsCMatrix.R, will give
130      faster determinant
131    
132    - tr(A %*% B) {and even  tr(A %*% B %*% C) ...} are also needed
133      frequently in some computations {conditional normal distr. ...}.
134      Since this can be done faster than by
135        sum(diag(A %*% B))  even for traditional matrices, e.g.
136                   sum(A * t(B)) or {even faster for "full" mat}
137                   crossprod(as.vector(A), as.vector(B))
138      and even more so for, e.g.  <sparse> %*% <dense>
139      {used in Soeren's 'gR' computations},
140      we should also provide a generic and methods.
141    
142    - qr.R(qr(x)) may differ for the "same" matrix, depending on it being
143      sparse or dense:
144        "qr.R(<sparse>) may differ from qr.R(<dense>) because of permutations"
145    
146      This is not really acceptable and currently influences  rcond() as well.
147    
148    - eigen() should become generic, and get a method at least for diagonal,
149      but also for symmetric -> dsyMatrix  [LAPACK dsyev() uses UPLO !],
150      but also simply for dgeMatrix (without going via tradition matrices).
151      What about Sparse?  There's fill-in, but it may still be sensible, e.g.
152      mlist <- list(1, 2:3, diag(x=5:3), 27, cbind(1,3:6), 100:101)
153      ee <- eigen(tcrossprod(bdiag(lapply(mlist, as.matrix))))
154      Matrix( signif(ee$vectors, 3) )
155    
156    - facmul() has no single method defined;  it looks like a good idea though
157      (instead of the infamous qr.qy, qr.qty,.... functions)
158    
159    - symmpart() and skewpart()  for *sparse* matrices still use (x +/- t(x))/2
160      and could be made more efficient.
161      Consider going via  asTuniq() or something very close to
162      .Arith.Csparse() in R/Ops.R
163    
164    - many setAs(*, "[dl]..Matrix") are still needed, as long as e.g.
165      replCmat() uses as_CspClass() and drop0(.) which itself call
166      as_CspClass() quite a bit.  --> try to replace these by
167      as(*, "CsparseMatrix"); forceSymmetric, etc.
168    
169    - implement fast diag(<triangularCsparse>) via calling new
170      src/Csparse.c's diag_tC_ptr()
171    
172    - add examples (and tests!) for update(<CHMfactor>, ..) and
173      Cholesky(......, Imult), also tests for hidden {hence no examples}
174      ldetL2up() { R/CHMfactor.R }
175    
176    - chol(<nsCMatrix>)  gives "temporarily disabled"
177      but should give the *symbolic* factorization;
178      similarly Cholesky(.) is not enabled
179    
180    - writeMM(obj, file=stdout()) creates file "1" since file is silently
181      assumed to be a string, i.e. cannot be a connection.
182      An R (instead of C) version should be pretty simple, and would work with
183      connections automatically ["lsparse" become either "real" or
184      "pattern", "depending if they have NAs or not].

Legend:
Removed from v.1087  
changed lines
  Added in v.2175

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge