search

scipy / scipy

SciPy library main repository
https://scipy.org
BSD 3-Clause "New" or "Revised" License
13.08k stars 5.19k forks source link

Add LAPACK's ZUNCSD cosine sine decomposition #11429

Closed balopat closed 4 years ago

balopat commented 4 years ago

Is your feature request related to a problem? Please describe.

I would like to use Cosine Sine decomposition in certain quantum compilation algorithms to decompose 3 qubit circuits - if this project of mine succeeds this will be implemented in Google's Cirq. CSD is Central to the algorithm. I do have a prototype working with https://github.com/anuesseler/pycsd however I believe SciPy would be a great place for this function.

Describe the solution you'd like

As scipy already has a mature process of importing and distributing compiled LAPACK functions, I would like to have ZUNCSD exposed through scipy.linalg. (Eventually maybe it would be great to have a higher level implementation of CSD amongst the other matrix decompositions.) I believe other python based quantum frameworks could benefit from having zuncsd exposed this way.

Describe alternatives you've considered

Alternatives considered:

Additional context (e.g. screenshots)

I am open to raise a PR for this (guidance is welcome if you are open to adding this).

ilayn commented 4 years ago

This sounds nice. We need to add also sorcsd, dorcsd, and cunscdas a group. @mdhaber Is this along with your lapack project?

balopat commented 4 years ago

@mdhaber, @ilayn, if you're okay to provide feedback I'm happy to give this a try and open a PR

mdhaber commented 4 years ago

No, this is not on the docket. Go for it! Thanks @balopat!

mdhaber commented 4 years ago

@balopat where did your most recent message go? I got an email but I don't see it here. If they need a separate wrapper and would be a fair amount more work, then it sounds OK to me for them to be in a separate PR.

balopat commented 4 years ago

Ah sorry, yes, I removed the comment as I found a bug in my test code. Maybe I should have just hidden it/edited it, sorry for the confusion! In the end it's not that bad, I can do them together. I have a working prototype for all four. Now I'm working on making the wrappers and the way I call the methods "scipy"-esque.

I have a couple of issues that I'd love to get guidance on:

Thank you!

ilayn commented 4 years ago

How would you suggest going about using multiple workspaces?

If you set all work parameters to -1 and expose them then you get all at once see for example ?heevd_lwork

https://github.com/scipy/scipy/blob/master/scipy/linalg/flapack_sym_herm.pyf.src#L195-L219

s/dor/c/zuncsd has IWORK - do we need to return that, or is it okay to intent(hide) it?

All work vars are hidden except for xxxx_lwork ones.

can you point me to a guide around this?

It is indeed a black magic but if you make an example we can have a look right away. You can basically follow the existing examples. Scalars and internal arrays almost always pass by reference.

balopat commented 4 years ago

heevd_lwork returns three arrays, however only one is used to pass in during the actual heevd call, lwork, the others are by default calculated?

https://github.com/scipy/scipy/blob/6bb954d3ce5c743b008bc9036dc29e6f34640b80/scipy/linalg/flapack_sym_herm.pyf.src#L182-L190

and here is the python side logic with _compute_lwork only returning lwork and only that is used in the lwork_args:

https://github.com/scipy/scipy/blob/6bb954d3ce5c743b008bc9036dc29e6f34640b80/scipy/linalg/decomp.py#L555-L560

Maybe I'm misinterpreting it, but it looks like heevd can "get away with this" because of the default values for lrwork 1+5*n+2*n*n:n and iwork. However, c/zuncsd don't have good defaults values for lrwork.

balopat commented 4 years ago

heevd_lwork returns three arrays, however only one is used to pass in during the actual heevd call, lwork, the others are by default calculated?

https://github.com/scipy/scipy/blob/6bb954d3ce5c743b008bc9036dc29e6f34640b80/scipy/linalg/flapack_sym_herm.pyf.src#L182-L190

and here is the python side logic with _compute_lwork only returning lwork and only that is used in the lwork_args:

https://github.com/scipy/scipy/blob/6bb954d3ce5c743b008bc9036dc29e6f34640b80/scipy/linalg/decomp.py#L555-L560

Maybe I'm misinterpreting it, but it looks like heevd can get away with this because of the default values for lrwork (1+5*n+2*n*n:n) and liwork (3+5*n:1). However, c/zuncsd don't have defaults at all for lrwork.

balopat commented 4 years ago

All work vars are hidden except for xxxx_lwork ones.

Ah, I just reread this, and now I see that both iwork/liwork and rwork/lrwork are hidden in heevd, that's how they "get away" with just passing in lwork, not the deafults. So, just for my curiosity - LWORK contains all the information necessary then for all the other workspaces?

ilayn commented 4 years ago

Maybe heevd was a bad example because those particular arrays have default and fixed sizes. Sorry about that. I'll find a better example where all of them are exposed in the lwork query and all hidden in the actual call.

ilayn commented 4 years ago

Let's try with evr for example

https://github.com/scipy/scipy/blob/6bb954d3ce5c743b008bc9036dc29e6f34640b80/scipy/linalg/flapack_sym_herm.pyf.src#L765-L772

Here you see that liwork and lwork are optional input vars, actual arrays are hidden. But if you look at its lwork counterpart it is the other way around

https://github.com/scipy/scipy/blob/6bb954d3ce5c743b008bc9036dc29e6f34640b80/scipy/linalg/flapack_sym_herm.pyf.src#L805-L810

balopat commented 4 years ago

Thanks @ilayn - I opened a Draft PR for some initial feedback, this is how I used _lwork:

https://github.com/scipy/scipy/pull/11524/files#diff-0547d172b04a54a4327aa6280727306cR130-R152

I'll start working my way through the TODO list, but any input is welcome!