[BUG/Feature Request] Matrix algebra

[fpavone]

It seems that RevBayes doesn't implement matrix algebra such as sum or product of two matrices:

A <- diagonalMatrix(2)
B <- diagonalMatrix(2)
A+B
   Error:   No overloaded function '_add' matches for arguments ( MatrixRealSymmetric, MatrixRealSymmetric )
    Potentially matching functions are:
       add (Natural<any> first,
            Natural<any> second)
       add (Integer<any> first,
            Integer<any> second)
       add (Real<any> first,
            Real<any> second)
       add (RealPos<any> first,
            RealPos<any> second)
       add (String<any> first,
            String<any> second)
       add (String<any> first,
            Real<any> second)
       add (String<any> first,
            Integer<any> second)
       add (String<any> first,
            Natural[]<any> second)
       add (String<any> first,
            Integer[]<any> second)
       add (String<any> first,
            RealPos[]<any> second)
       add (String<any> first,
            Real[]<any> second)
       add (Natural[]<any> first,
            Natural[]<any> second)
       add (Integer[]<any> first,
            Integer[]<any> second)
       add (RealPos[]<any> first,
            RealPos[]<any> second)
       add (Real[]<any> first,
            Real[]<any> second)
       add (Natural<any> first,
            Natural[]<any> second)
       add (Integer<any> first,
            Integer[]<any> second)
       add (Real<any> first,
            Real[]<any> second)
       add (RealPos<any> first,
            RealPos[]<any> second)

I think it could be useful to include matrix sum and product in RevBayes.

Thanks

[fpavone]

Actually, it looks that internally this operations are defined in the header related to MatrixReal src/core/dataypes/math/MatrixReal.h, but it seems they are not "exposed" to the language. Is it the case or I'm missing something?

Moreover, it seems that it's not even possible to define a user-constructed MatrixReal. For example, supposed that I want to create a random variable with a multivariate normal distribution with a specific covariance matrix. How can I define a MatrixRealSymmetric to specify my covariance matrix?

I struggled to find any hint about this from the documentations and the tutorials.

[bredelings]

It looks like none of the developers have needed this themselves, and so they haven't added it.

Normally, to add a function at the Rev language level, you need to

1. create a Rev function (for example src/revlanguage/functions/math/Func_sqrt.cc)
2. add it to src/revlanguage/workspace/RbRegister_Func.cc

However, in this case it seems like there are already templates to make the C++ functions available at the Rev language level in RbRegister_Basic.cpp. I will see if I can submit a PR for this.

I suspect that we might need to add a way to create a constant matrix. It looks like people have mostly been drawing matrices from priors, and therefore haven't added a way to create a constant matrix. I will ask on the slack group...

[bredelings]

Hmm... adding the matrix operations is going to be a bit more tricky than I thought.

[bredelings]

OK, I figured out how to add addition, subtraction, and negation.

[bredelings]

You can make your own matrices like so:

> sigma = matrix( v( v(1,0), v(0,1) ) )
> sigma
   [ [ 1.0000, 0.0000 ] ,
     0.0000, 1.0000 ] ]
> type(sigma)
   MatrixReal

However, these matrices have type MatrixReal, and not MatrixRealSymmetric. So they can't be used as covariance matrices (yet).

[fpavone]

Thank you very much for the fast reply!

OK, I figured out how to add addition, subtraction, and negation.

This is great! I think having some basic matrix algebra (as also matrix multiplication) available can make RevBayes even more flexible.

You can make your own matrices like so:

sigma = matrix( v( v(1,0), v(0,1) ) ) sigma [ [ 1.0000, 0.0000 ] , 0.0000, 1.0000 ] ] type(sigma) MatrixReal

However, these matrices have type MatrixReal, and not MatrixRealSymmetric. So they can't be used as covariance matrices (yet).

My bad I didn't see it in the documentation. I hope there will be soon also the possibility to create MatrixRealSymmetric (or cast it from MatrixReal), making it possible to use user-made matrices as input to functions requiring specific types. Maybe the check for symmetric matrices could be done internally to functions requiring it, so to have these functions requiring just a MatrixReal in input. Although, I understand that it may be non-trivial to do it properly. Also I guess memory efficiency can be a reason of why having a specific type for symmetric matrices.