Re: [stan-dev/math] Create basic lower-triangular matrix type (#364)

45 views
Skip to first unread message

Bob Carpenter

unread,
Sep 2, 2016, 7:50:19 AM9/2/16
to stan...@googlegroups.com
[off GitHub]

I'm trying to keep our GitHub issue discussions minimal and focused
on concrete deliverables.

Adding a new type is a substantial amount of work in the stan-dev/math
library (all the functions and everything to manipulate type
and transforms), in stan-dev/stan (all the parsing and code generation
as well as all of the services and I/O, which includes all the type
inference rules), and in all of the interfaces (implementing the generic
callbacks from the services to handle I/O w.r.t. Stan and then all of the
external output).


DESIGN

Sounds like you're committed to adding a new basic type
(alongside int, real, vector, row_vector, and matrix). That's
the long road to getting this done and is going to be confusing
to users, but is going to be necessary to squeeze every bit
of performance out of it. Then you have to decide how it's going
to interact with regular matrices, and in particular if mappings
to matrix types are implicit or explicit (matrix to tri_matrix would
have to be explicit, but the other way around can be implicit).


MATH LIB

You need to decide on the underlying Eigen type you want to
use to represent the lower-triangular matrix type. I think
the decision is whether you want to use the view template expression
or an underlying dense matrix or just an array/vector to store
the values. Unless the answer is dense matrix with zeros above
diagonal, you'll have to make sure all the functions to which you
want it to apply to accept your type. If you do use the dense
formulation (which may make your life easier downstream), make sure
all the zero entries are copies of the same underlying var or fvar.

Big question is whether the number of rows and columns must match
or not.

You'll need the transforms for constraining and unconstraining (even
if it's just the identity map) with log Jacobians (may be zero).


LANGUAGE

Basically, you need to deal with the parser and code generator
and plumb through the math lib

Parser

You're right that you need to declare a new variable declaration
parser rule with semantics (_def.hpp) and declaration of the semantic
type (.hpp). You'll need a new semantic type (like DOUBLE_T).

But that's just the tip of the iceberg. The much more involved
part is all the type inference inside the language that computes
the type and dimensions of expressions. And it figures out what
can be done for assignment and function calls. For example, can
you send in a lower triangular matrix into a regular matrix-accepting
function? You can pass an int type into a real function, but that's
the only built-in promotion we have. You'll see lots of ad hoc logic
throughout the parser semantic actions for dealing with this.

Code Generator

You need to make sure the efficient version of the math lib types
get generated where they're needed. How hard this is will depend on
the choices in the math library and for conversion.


I/O

You'll need to figure out how to use var_context to read and
how the writer callbacks will write.


SERVICES

Not sure what you need to do here to get the I/O to the
services. Maybe nothing.


INTERFACES

Then you need to plumb it through the interfaces. For example,
CmdStan will need a dump format (if we keep that format, then
the services refactor won't matter) and RStan and PyStan will need
to figure out how to write all the appropriate callbacks.



That's all off the top of my head. There are going to be many
more little details that come up for all of this.


- Bob


> On Sep 2, 2016, at 2:16 AM, Rob Trangucci <notifi...@github.com> wrote:
>
> @bob-carpenter I talked to @syclik yesterday about whether any of the refactor would interfere with the work needed to pipe a new type through the language and he didn't think it would. I'd like to start working on it this weekend; do you have an idea of what changes I'd need to make to src/stan/lang/generator.hpp or src/stan/lang/ast.hpp to get a new matrix type working? Right now it seems like I'll need to:
>
> • write a lower_triangular_var_decl similar to cholesky_corr_var_decl
> • modify write_array_vars_visgen to handle lower_triangular_var_decl
> I'm sure there's more, but I'm not sure where to go from here.
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub, or mute the thread.
>

Ben Goodrich

unread,
Sep 2, 2016, 10:14:49 AM9/2/16
to stan development mailing list
On Friday, September 2, 2016 at 7:50:19 AM UTC-4, Bob Carpenter wrote:
You need to decide on the underlying Eigen type you want to
use to represent the lower-triangular matrix type.  I think
the decision is whether you want to use the view template expression
or an underlying dense matrix or just an array/vector to store
the values.  

Also whether to store the vector that represents the lower triangular part of the matrix in row-major or column-major order. For most of the operations people tend to do with lower triangular matrices, I think row-major is going to be more computationally efficient.

Ben
Reply all
Reply to author
Forward
0 new messages