matrix logarithm function for mxAlgebra

9 replies [Last post]
CharlesD's picture
Offline
Joined: 04/30/2013

Any chance of this being implemented at some point? Would seem to be consistent with the omxExponential implementation. Currently I optimize some free parameters over the range 0 to -inf, which is exponentially scaled, and I suspect this probably poses optimization difficulties....

jpritikin's picture
Offline
Joined: 05/23/2012
matrix logarithm

If you can point me to an open-source C or C++ implementation that we can include, I'll add it in a jiffy.

mhunter's picture
Offline
Joined: 07/31/2009
Confusion

There might be some confusion. We currently have a log() function for mxAlgebras. Inside an mxAlgebra, log(aMatrixOrAlgebra) will produce the elementwise logarithm of aMatrixOrAlgebra. The same is true for exp(aMatrixOrAlgebra): it's an elementwise exponential of aMatrixOrAlgebra. The omxExponential() function is an entirely different thing. It is the matrix exponential function, not the elementwise exponential of a matrix. Similarly, the matrix logarithm function is not the elementwise logarithm function.

If you want to take the log of every element in an algebra or matrix, you can do that with log() already. If you want to undo the matrix exponential, you'll need a matrix logarithm. Two issues with matrix logarithms in general are that (1) not all square matrices have a logarithm, and (2) among those that do, some have more than one.

Your question sounded like you wanted the elementwise logarithm that we already have. Or did you actually want the matrix logarithm?

jpritikin's picture
Offline
Joined: 05/23/2012
matrix log

Charles is the guy working on continuous time models. I'm pretty sure he want the matrix log, not the elementwise log.

I assume there is some way to constrain the matrix log such that there is only 1 solution.

tbates's picture
Offline
Joined: 07/31/2009
some sources
CharlesD's picture
Offline
Joined: 04/30/2013
Joshua is correct, it's the

Joshua is correct, it's the matrix log I'd like. I'm not clear under which conditions it is not defined, (something I should get clear on, certainly!) but I ran a limited range of tests and logm() from the expm package gave me consistent results. If it proves inconsistent then I could just leave it as optional...

A side issue - does my intuition that the exponential scale of the free parameters may be posing estimation difficulties in some circumstances seem reasonable?

jpritikin's picture
Offline
Joined: 05/23/2012
optimizer

> A side issue - does my intuition that the exponential scale of the free parameters may be posing estimation difficulties in some circumstances seem reasonable?

I think so. If the true value of a parameter is less than about 1e-9 then the optimizer will believe it is done because changes to that value are already smaller than the tolerance. Log units don't have this problem.

jpritikin's picture
Offline
Joined: 05/23/2012
try it out

SVN 3592

Let us know how it goes.

I used the algorithm with a C implementation. There is a better algorithm available, Higham 2008, but we need to translate it from R to C.

CharlesD's picture
Offline
Joined: 04/30/2013
Works well, though yes there

Works well, though yes there are some more complex cases I deal with where the higham08 alg would cope but the simple approach doesn't. Need to properly test to see if optimization for me has improved in general, but cheers! :)

jpritikin's picture
Offline
Joined: 05/23/2012
higham08

It would probably take me 1-2 days to translate the R code into C. It's not a huge amount of work, but there are lots of other demands on my time. Let me know if you think it would really be worth it.