-
Notifications
You must be signed in to change notification settings - Fork 2.5k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Op/optimization TODO: log determinant #150
Comments
Related: there is an optimization in linalg for the log det of PSD On Mon, Oct 24, 2011 at 11:46 AM, David Warde-Farley
|
David, can you be more explicit in what is needed? A new op that call numpy.slogdet and an optimization that replace log(abs(det(X)) by it? We need to do some speed comparison before enabling it by default. |
This would be a stability optimization, the speed shouldn't prevent On Mon, Oct 24, 2011 at 11:57 AM, nouiz
|
I would be interested to know how much slower it is. Having a stability optimization that is 10 times slower is not the best. In that case, we should document it well and tell people that don't want to how to disable it. |
We should be able to replace log(det(X)) with an slogdet Op as well, and just raise an error if the sign is negative (I guess this would be accomplished by an |
@nouiz: Newer versions of NumPy actually call We have an old version of NumPy at the lab which does not include I actually just did a comparison of the old In [28]: timeit linalg.slogdet(X)
100 loops, best of 3: 15.1 ms per loop
In [29]: timeit det(X) # old implementation, copied from NumPy 1.4.1
100 loops, best of 3: 15.4 ms per loop |
@jaberg That raises the interesting issue of how to know that certain intermediate results are positive-definite matrices, and how to mark them as such. Seems like something that the user would need to specify in most cases, like Any thoughts on the right way to do this? |
So for the speed, it seam good. James implemented an Hint Env feature to allow the user to give hint to the compiler. He used it for the other optimization he talked about. This approch seam good for me, but it can take much time to implement correctly all case correctly. At least that is what took time with the ShapeFeature. |
On Mon, Oct 24, 2011 at 2:29 PM, David Warde-Farley
Check out how it's done in sandbox.linalg. there is a function
|
Just out of curiosity, why hasn't this been done yet? I haven't coded any Theano ops, but I could give it a shot... maybe just a log abs det, since that's probably the most common case (in MLE estimation etc)? FYI I'm hitting the stability issue pretty badly... |
I'm not sure why this has not been done yet, or maybe some people ended up implementing that in their own repositories. |
OK I'll give it a try! On Wed, Jan 27, 2016, 22:39 Pascal Lamblin notifications@github.com wrote:
|
Found a couple of implementations on Github. This one does an SVD first and takes the log of the diagonal (not sure why it's diagonal squared... need to check the math; also has a funny expression for the gradient... it's correct but computationally expensive). This one uses numpy's slogdet. Any preferences? SVD on GPU should be a simple copy/paste job I think, but I haven't found any CUDA/PyCUDA/gnumpy code for the slogdet... |
You could check inside numpy code how slogdet is implemented. Maybe it does There is svd code on the GPU somewhere, so with that, it wouldn't be too On Thu, Jan 28, 2016 at 3:25 AM, Heikki Arponen notifications@github.com
|
It calls something in lapack routines in numpy/numpy/linalg/umath_linalg.c.src, which says in docstring that it "computes logdet from factored diagonal" so I guess it is SVD (or really an eigendecomposition since it's for square matrices)... The square in the svd code was wrong as I suspected, but fixing that, the SVD and slogdet versions give slightly different results! Difference of about 1E-5 for 100*100 random normal matrices about 80% of time and zero difference 20% of time... happens for pure numpy/scipy versions too. Maybe there's some numerical stability issues in first doing the svd, then taking log and then sum?? I'll post the code below if anyone wants to take a look. I'm actually taking log abs det, because why take a log of a negative number... So I'm guessing the slogdet version is numerically more accurate, which raises the question of how to do the SVD version accurately on a GPU... dammit! PS. This will (hopefully) be my first ever PR, so I need some time to figure out how to code the tests, where to put them etc... Numpy slogdet:
The SVD is the same, except now there's this inside the perform:
|
For float64 the difference is 1E-14:
For float32 the difference is 1E-5!
Anyone know the reason for this? EDIT: craaaap something wrong with my numpy/scipy installation... can anyone test if you get the same errors? Probably not... need to reinstall/recompile numpy and scipy I guess... |
Given that the logdet values for the examples above are of order 1E+2, the relative error is around 1E-7...1E-8 for float32... maybe that's acceptable? I think I should go with the direct SVD method then for the 'perform' method, if the GPU version is going to be also SVD... |
Yes, this kind of differences is really in the normal range for float32, these are even quite low. |
OK did a PR here: #3954 By the way, SVD on GPU seems to be useful only when matrix dims go above 1000 * 1000, so maybe the GPU implementation is not that important yet? Or will there be some overhead from loading parameters back and forth between GPU memory and RAM?? This image is from this paper (not very fresh though...): |
New PR due to me being a n00b: #3959 |
The transfer to/from the GPU/CPU is to be avoided as much as possible. So On Sun, Jan 31, 2016 at 3:45 PM, Heikki Arponen notifications@github.com
|
OK got it |
Now that we have a determinant Op we should probably look into this.
There are numerically stable ways to compute the natural logarithm of the determinant, a quantity which is often needed in evaluating the log probability of a multivariate Gaussian distribution. NumPy provides this as
numpy.slogdet
for "sign log determinant" (it computeslog(abs(det(X))
) and also returns the sign of the determinant).The derivative of the log determinant also has a particularly simple form as
inv(X).T == inv(X.T)
.The text was updated successfully, but these errors were encountered: