[CppAD] Format for sparse Jacobians/Hessians

[Michael Braun]

Bob:

I apologize for not being more specific in my proposal.  Yes, sparse_hessian would be a member function of the ADFun<Base> object.  I really envisioned two functions:  one that would return the sparsity structure, and another that would take that structure as input, and return the values.  One might already know the structure, and certainly one might want to save the structure from previous calls for the same matrix.  And of course, there could be a single function that wraps together computing the structure and the values.

As far as the output format goes, I had not considered which format might be easiest to compute.  If triplets are the easiest, then that's fine with me, as long as the output is sorted, and that the user has an option for row or column ordering, and lower or upper triangle.  Converting from triplets to CSC, or many other formats, is straightforward, although it would involve some additional overhead.

The "big picture" objective is to avoid allocating or returning an NxN dense matrix, even internally, when computing sparse Hessians, which is what CppAD does now.

Michael



> 
> Michael et al.,
> 
> Your proposed get_sparse_hessian interface seems backward to me. CppAD
> can figure out which elements are nonzero; the application code calling
> it can't, unless it does CppAD's work for it. Hence the indices defining
> the sparse structure of the matrix should be outputs, not inputs (not
> const). Am I missing something?
> 
> Apparently get_sparse_hessian is to be a member function of class
> ADFun<Base>. If not, it needs another argument, the ADFun object to be
> evaluated.
> 
> 
> describes a large number of sparse matrix formats that are efficient for
> input to sparse solvers and other standard operations. The trouble with
> all of them is that the function that fills them must proceed in strict
> row-major or column-major order. Backtracking would be very expensive.
> Now, I don't know exactly how the innards of a function like
> get_sparse_hessian would work within CppAD, but I suspect that elements
> would get built up as they are encountered in the tape, not in any
> particular index order. If that's true, it seems to me that the
> evaluation pretty much has to be in terms of triplets. I have used a
> std::map<std::pair<size_t,size_t>,Value_type> successfully, with
> reasonable speed.
> 
> Admittedly, a "map" format is very poor for most subsequent matrix
> operations. That is often the case: Structures that are good for one
> phase of a data-set's lifetime are poor for another phase. But that's
> OK, because traversing a std::map in row or column order is easy and
> fast, so it can be transformed into any of the "standard" formats trivially.
> 
> The focus of CppAD is the first phase, creating the sparse matrix. The
> best choice of the representation for the next phase depends critically
> on what the client software needs to do next. Hence, unless that next
> operation is to be part of the CppAD library, CppAD should stay out of
> that debate. Choose a "filling" format that is efficient for filling the
> matrix, and let the application code transform it to a "working" format
> if that makes sense, in preparation for whatever comes next.






-------------------------------------------
Michael Braun
Associate Professor of Management Science (Marketing Group)
MIT Sloan School of Management
100 Main St.., E62-535
Cambridge, MA 02139
braunm at mit.edu
617-253-3436




-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://list.coin-or.org/pipermail/cppad/attachments/20120110/3cc6343e/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/pkcs7-signature
Size: 1844 bytes
Desc: not available
URL: <http://list.coin-or.org/pipermail/cppad/attachments/20120110/3cc6343e/attachment.p7s>