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

LogAbsDet #3959

Open
wants to merge 6 commits into
from

Conversation

Projects
None yet
7 participants
@harpone

harpone commented Jan 31, 2016

Travis CI corrections (Python 3.3 compatibility fixes; PEP8 fixes)

Second go at the LogAbsDet op. Has gradient, but no GPU code so far... if anyone wants to help/ point me to PyCuda/CUDA implementation of SVD, let me know!

This is very similar to numpy.linalg.slogdet, except that the sign of the determinant is omitted, so it just computes log(abs(det(M))) for a matrix M.

Note that M does not need to be positive semidefinite, since we're taking the absolute value (as happens in e.g. maximum likelihood estimation).

do the op for gh-150, the optimization isn't done.

(x,) = inputs
(z,) = outputs
s = numpy.linalg.svd(x, compute_uv=False)
log_abs_det = numpy.sum(numpy.log(numpy.abs(s)))

This comment has been minimized.

@nouiz

nouiz Jan 31, 2016

Member

I didn't check the math. But this can be done by reusing current op. The svd op and the elemwise op. This will be better as the elemwise op will get optimized by Theano and will work on the GPU (for the elemwise).

Then the only think missing to make it work on the GPU would be to have an svn op on the GPU. There is one implementation here that won't request c code: http://scikit-cuda.readthedocs.org/en/latest/generated/skcuda.linalg.svd.html

@nouiz

nouiz Jan 31, 2016

Member

I didn't check the math. But this can be done by reusing current op. The svd op and the elemwise op. This will be better as the elemwise op will get optimized by Theano and will work on the GPU (for the elemwise).

Then the only think missing to make it work on the GPU would be to have an svn op on the GPU. There is one implementation here that won't request c code: http://scikit-cuda.readthedocs.org/en/latest/generated/skcuda.linalg.svd.html

This comment has been minimized.

@nouiz

nouiz Feb 1, 2016

Member

In fact, I think it is a good idea to keep this as an op. It will lower memory usage during the gradient on the CPU and will provide a gradient (the currrent svd op don't have a grad implemented).

but an easy way to have this on the GPU would be to make an optimization that will convert this op to a graph that reuse the elemwise on the GPU and make an svd op. For the gradient, MatrixInverse will also need a GPU op to have this work on the GPU. But the first step, is to get this working!

@nouiz

nouiz Feb 1, 2016

Member

In fact, I think it is a good idea to keep this as an op. It will lower memory usage during the gradient on the CPU and will provide a gradient (the currrent svd op don't have a grad implemented).

but an easy way to have this on the GPU would be to make an optimization that will convert this op to a graph that reuse the elemwise on the GPU and make an svd op. For the gradient, MatrixInverse will also need a GPU op to have this work on the GPU. But the first step, is to get this working!

This comment has been minimized.

@harpone

harpone Feb 1, 2016

Yeah, I have no idea what the gradient of SVD would be :) grad of log abs det is very simple on the other hand, so maybe it's better to leave this out of the theano optimizations...

OK didn't know there's no GPU op for MatrixInverse... so step 1 would now be to get LogAbsDet work on GPU? Is it possible to use the skcuda.linalg.svd directly? (I have no idea how the make_thunk works...). So not sure if I should try to do the GPU SVD here directly or in the SVD op? Damn it's still a bit confusing :)

@harpone

harpone Feb 1, 2016

Yeah, I have no idea what the gradient of SVD would be :) grad of log abs det is very simple on the other hand, so maybe it's better to leave this out of the theano optimizations...

OK didn't know there's no GPU op for MatrixInverse... so step 1 would now be to get LogAbsDet work on GPU? Is it possible to use the skcuda.linalg.svd directly? (I have no idea how the make_thunk works...). So not sure if I should try to do the GPU SVD here directly or in the SVD op? Damn it's still a bit confusing :)

This comment has been minimized.

@nouiz

nouiz Feb 5, 2016

Member

We didn't document much making new gpu op as we are in transition to a new gpu back-end...

For a first pass, you can make one GPU op that have its gradient on the CPU and use it directly. Then after, you can add an optimization that will convert the CPU op to the GPU op (I can help, it is very easy when you know how to do it).

To make a GPU op, I would use as an example the code of CuFFTOp in sandbox/cuda/fftconv.py. This use scikits.cuda for computation, so you can reuse it for the svd implementation.

@nouiz

nouiz Feb 5, 2016

Member

We didn't document much making new gpu op as we are in transition to a new gpu back-end...

For a first pass, you can make one GPU op that have its gradient on the CPU and use it directly. Then after, you can add an optimization that will convert the CPU op to the GPU op (I can help, it is very easy when you know how to do it).

To make a GPU op, I would use as an example the code of CuFFTOp in sandbox/cuda/fftconv.py. This use scikits.cuda for computation, so you can reuse it for the svd implementation.

This comment has been minimized.

@harpone

harpone Feb 5, 2016

OK thanks, I'll probably get that done some time soon (been a bit busy lately)!

@harpone

harpone Feb 5, 2016

OK thanks, I'll probably get that done some time soon (been a bit busy lately)!

f = theano.function([x], self.op(x))
out = f(input_mat)
svd_diag = numpy.linalg.svd(input_mat, compute_uv=False)
numpy_out = numpy.sum(numpy.log(numpy.abs(svd_diag)))

This comment has been minimized.

@nouiz

nouiz Jan 31, 2016

Member

I think it would be better to compare again the numpy function you want to implement.

@nouiz

nouiz Jan 31, 2016

Member

I think it would be better to compare again the numpy function you want to implement.

This comment has been minimized.

@harpone

harpone Feb 1, 2016

That would be the slogdet function then, but it probably does something a bit different since the results are differ (relatively) by about 1E-8, as I discussed in the original issue (here)[https://github.com/Theano/Theano/issues/150]. Maybe I should compare against the slogdet but increase the tolerance of the allclose?

@harpone

harpone Feb 1, 2016

That would be the slogdet function then, but it probably does something a bit different since the results are differ (relatively) by about 1E-8, as I discussed in the original issue (here)[https://github.com/Theano/Theano/issues/150]. Maybe I should compare against the slogdet but increase the tolerance of the allclose?

This comment has been minimized.

@nouiz

nouiz Feb 5, 2016

Member

Having 1e-8 difference isn't a problem. You can have that type of error just by changing the order of addition when you add a vector of floats.

If the utt.assert_allclose pass, the I would compare again numpy.slogdet. If not, we need to understand why as they are supported to be equivalent. Maybe the fix would be to raise the comparison tolerance.

@nouiz

nouiz Feb 5, 2016

Member

Having 1e-8 difference isn't a problem. You can have that type of error just by changing the order of addition when you add a vector of floats.

If the utt.assert_allclose pass, the I would compare again numpy.slogdet. If not, we need to understand why as they are supported to be equivalent. Maybe the fix would be to raise the comparison tolerance.

This comment has been minimized.

@harpone

harpone Feb 5, 2016

I think my test cases are a bit silly, as the svd values are pretty big, about 1e+2 or so.. so maybe I should just rescale the test matrices (by e.g. 1 / sqrt(size)), then the error will probably be below the tolerance.

@harpone

harpone Feb 5, 2016

I think my test cases are a bit silly, as the svd values are pretty big, about 1e+2 or so.. so maybe I should just rescale the test matrices (by e.g. 1 / sqrt(size)), then the error will probably be below the tolerance.

This comment has been minimized.

@nouiz

nouiz Feb 24, 2016

Member

Did you looked into this? I'm checking if we can merge this PR.

@nouiz

nouiz Feb 24, 2016

Member

Did you looked into this? I'm checking if we can merge this PR.

This comment has been minimized.

@harpone

harpone Feb 24, 2016

Sorry, haven't had the time to look into it as I just started a new day job recently... but it's definitely on my to-do list, just need to find a window for it

@harpone

harpone Feb 24, 2016

Sorry, haven't had the time to look into it as I just started a new day job recently... but it's definitely on my to-do list, just need to find a window for it

@lamblin

This comment has been minimized.

Show comment
Hide comment
@lamblin

lamblin Feb 4, 2016

Member

Forcing new run of Travis after #3885 has been merged.

Member

lamblin commented Feb 4, 2016

Forcing new run of Travis after #3885 has been merged.

@lamblin lamblin closed this Feb 4, 2016

@lamblin lamblin reopened this Feb 4, 2016

@twiecki

This comment has been minimized.

Show comment
Hide comment
@twiecki

twiecki Mar 5, 2016

This might be a relevant PR for pymc3 too. Any way that the theano simplification would change a manual log(det()) to use this logdet directly?

twiecki commented Mar 5, 2016

This might be a relevant PR for pymc3 too. Any way that the theano simplification would change a manual log(det()) to use this logdet directly?

Show outdated Hide outdated theano/tensor/nlinalg.py
raise
def grad(self, inputs, g_outputs):
gz, = g_outputs

This comment has been minimized.

@superbobry

superbobry Mar 5, 2016

Contributor

The comma here can be hard to spot, why not use single-element list syntax instead?

[gz] = g_outputs
@superbobry

superbobry Mar 5, 2016

Contributor

The comma here can be hard to spot, why not use single-element list syntax instead?

[gz] = g_outputs

This comment has been minimized.

@harpone

harpone Mar 5, 2016

[gz] = g_outputs

yeah sounds good...

@harpone

harpone Mar 5, 2016

[gz] = g_outputs

yeah sounds good...

This comment has been minimized.

@harpone

harpone Mar 5, 2016

Any way that the theano simplification would change a manual log(det()) to use this logdet directly?

There may be some problems with that if positivity of det is not enforced... not sure how to handle that...

@harpone

harpone Mar 5, 2016

Any way that the theano simplification would change a manual log(det()) to use this logdet directly?

There may be some problems with that if positivity of det is not enforced... not sure how to handle that...

@superbobry

This comment has been minimized.

Show comment
Hide comment
@superbobry

superbobry Mar 5, 2016

Contributor

This looks great! Is it possible to automatically optimize log(det(...)) to logdet(...)?

Contributor

superbobry commented Mar 5, 2016

This looks great! Is it possible to automatically optimize log(det(...)) to logdet(...)?

s = numpy.linalg.svd(x, compute_uv=False)
log_abs_det = numpy.sum(numpy.log(numpy.abs(s)))
z[0] = numpy.asarray(log_abs_det, dtype=x.dtype)
except Exception:

This comment has been minimized.

@superbobry

superbobry Mar 5, 2016

Contributor

Is it possible to be more specific about the exception type here?

@superbobry

superbobry Mar 5, 2016

Contributor

Is it possible to be more specific about the exception type here?

This comment has been minimized.

@harpone

harpone Mar 5, 2016

Is it possible to be more specific about the exception type here?

hmm I guess one should look inside numpy.linalg.svd for similar clauses... or the tests??

Actually I'm pretty swamped with other stuff at the moment, so I'd appreciate if someone else finished this stuff :)

@harpone

harpone Mar 5, 2016

Is it possible to be more specific about the exception type here?

hmm I guess one should look inside numpy.linalg.svd for similar clauses... or the tests??

Actually I'm pretty swamped with other stuff at the moment, so I'd appreciate if someone else finished this stuff :)

This comment has been minimized.

@harpone

harpone Mar 5, 2016

This looks great! Is it possible to automatically optimize log(det(...)) to logdet(...)?

see here

@harpone

harpone Mar 5, 2016

This looks great! Is it possible to automatically optimize log(det(...)) to logdet(...)?

see here

Smart scaling of logabsdet tests
minor semantic fixes
@harpone

This comment has been minimized.

Show comment
Hide comment
@harpone

harpone Mar 5, 2016

Fixed the tests and did the minor change as suggested by @superbobry but Travis CI seems to be failing... don't know why and I have no time to check now

harpone commented Mar 5, 2016

Fixed the tests and did the minor change as suggested by @superbobry but Travis CI seems to be failing... don't know why and I have no time to check now

@abergeron abergeron closed this Mar 7, 2016

@abergeron abergeron reopened this Mar 7, 2016

@abergeron

This comment has been minimized.

Show comment
Hide comment
@abergeron

abergeron Mar 7, 2016

Member

There was some breakage in the travis setup. A the new run should go fine (assuming there are no problems in this PR).

Member

abergeron commented Mar 7, 2016

There was some breakage in the travis setup. A the new run should go fine (assuming there are no problems in this PR).

@harpone

This comment has been minimized.

Show comment
Hide comment
@harpone

harpone Mar 12, 2016

@nouiz I think this is done now (except for the GPU version)

harpone commented Mar 12, 2016

@nouiz I think this is done now (except for the GPU version)

@lamblin

This comment has been minimized.

Show comment
Hide comment
@lamblin

lamblin Mar 23, 2016

Member

Now that the release is over, is there anything preventing this from being merged?

Member

lamblin commented Mar 23, 2016

Now that the release is over, is there anything preventing this from being merged?

@nouiz

This comment has been minimized.

Show comment
Hide comment
@nouiz

nouiz Mar 23, 2016

Member

I have comments that haven't been done. Mostly, the author started a new jobs, so he have difficulty to find time to finish this.

Member

nouiz commented Mar 23, 2016

I have comments that haven't been done. Mostly, the author started a new jobs, so he have difficulty to find time to finish this.

@harpone

This comment has been minimized.

Show comment
Hide comment
@harpone

harpone Mar 24, 2016

It should be working now. I fixed the tests too 12 days ago. The only thing missing is the CUDA code, which is of course a major omission... prolly can't find time for that any time soon...

harpone commented Mar 24, 2016

It should be working now. I fixed the tests too 12 days ago. The only thing missing is the CUDA code, which is of course a major omission... prolly can't find time for that any time soon...

@springcoil

This comment has been minimized.

Show comment
Hide comment
@springcoil

springcoil Jan 28, 2017

Is it possible to get this merged in without the CUDA code, and seperate that to another PR?
Due to it's interest for PyMC3

Is it possible to get this merged in without the CUDA code, and seperate that to another PR?
Due to it's interest for PyMC3

@twiecki

This comment has been minimized.

Show comment
Hide comment
@twiecki

twiecki Feb 16, 2017

@nouiz What would adding GPU support entail? Writing CUDA? Presumably SVD already exists for GPU.

twiecki commented Feb 16, 2017

@nouiz What would adding GPU support entail? Writing CUDA? Presumably SVD already exists for GPU.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment