PCA for short-fat arrays

[eric-czech]

Assuming I wanted to run non-randomized PCA on a short and fat array (chunks only along columns), is there a way to do that using the current dask-ml PCA implementation?

[TomAugspurger]

The memory usage is outlined in the notes at https://docs.dask.org/en/latest/array-api.html#dask.array.linalg.tsqr. So it should finish as long as the chunks are too big, it might just be slow?

[eric-czech]

Thanks @TomAugspurger, but this is what I mean (i.e. PCA.fit won't define the arrays):

from dask_ml.decomposition import PCA
# Chunk along columns but not rows
x = da.random.random(size=(100, 1000)).rechunk((-1, 50))
PCA(svd_solver='full').fit(x) # full or tsqr calls da.linalg.svd
...
/opt/conda/lib/python3.7/site-packages/dask_ml/decomposition/pca.py in _fit(self, X)
    272 
    273         if solver in {"full", "tsqr"}:
--> 274             U, S, V = da.linalg.svd(X)
    275         else:
    276             # randomized
...
ValueError: Input must have the following properties:
  1. Have two dimensions
  2. Have only one column of blocks

Note: This function (tsqr) supports QR decomposition in the case of
tall-and-skinny matrices (single column chunk/block; see qr)Current shape: (100, 1000),
Current chunksize: (100, 50)

I could call da.linalg.svd myself with the transposed matrix input and use whatever formulas there are that equate the singular vectors for a transposed input (iirc it's possible), but the current PCA code wouldn't work as is. I wanted to see if this has come up before though I don't see any related issues. It's a fairly common use case in genetics (cc: @mrocklin).

[mrocklin]

I could call da.linalg.svd myself with the transposed matrix input and use whatever formulas there are that equate the singular vectors for a transposed input (iirc it's possible), but the current PCA code wouldn't work as is

This would be my first instinct. (indeed, we might put this logic into da.linalg.svd directly?)

Are there problems that you ran into, or that you think you would run into with this approach?

[eric-czech]

Are there problems that you ran into, or that you think you would run into with this approach?

No not really except that some PCA implementation is available in the general domain for in memory arrays, tall-skinny ones, or fully chunked ones, but short-fat ones are just common enough (as will be the desire to avoid non-deterministic results) that I'd like to find a good solution for them. It will be awkward to write something into sgkit for just that without supporting all of them and I'd like to see us avoid it if we can.

In principle, I guess it's as simple as this:

import dask.array as da
import numpy as np 

x = da.random.random(size=(1000, 100)).rechunk((-1, -1))

# Run on tall-skinny array
u, s, v = da.linalg.svd(x)
pcs1 = v.compute()

# Run on transposition to short-fat array (USV becomes VtSUt)
vt, s, ut = da.linalg.svd(x.T)
pcs2 = vt.T.compute()

# Make sure the PCs are the same (modulo signs)
all([
    np.allclose(pcs1[i], pcs2[i]) or 
    np.allclose(pcs1[i], -pcs2[i]) 
    for i in range(x.shape[1])
])
> True

There are several other references to shapes and axes in the current PCA implementation that would keep this from working as-is, even if it is that simple to get the original components using a transposed input.

If there was a way for it to be handled transparently in da.linalg.svd, that would be terrific. It wasn't clear to me in a glance at that code why compute_svd is an option on tsqr but not sfqr. Maybe it would be trivial otherwise?

[mrocklin]

I have a hazy memory of trying this logic, but finding that singular values/vectors switched sign when I transposed, and so bailed on the effort due to not wanting to think about it. I could be fabricating that memory though. In any event, I'm +1 if anyone wants to pick this up. My guess is that it would be a small code contribution with some nice effects.

…

On Thu, Aug 27, 2020 at 2:58 PM Eric Czech ***@***.***> wrote: Are there problems that you ran into, or that you think you would run into with this approach? No not really except that some PCA implementation is available in the general domain for in memory arrays, tall-skinny ones, or fully chunked ones, but short-fat ones are just common enough (as will be the desire to avoid non-deterministic results) that I'd like to find a good solution for them. It will be awkward to write something into sgkit for just that without supporting all of them and I'd like to see us avoid it if we can. In principle, I guess it's as simple as this: import dask.array as daimport numpy as np x = da.random.random(size=(1000, 100)).rechunk((-1, -1)) # Run on tall-skinny arrayu, s, v = da.linalg.svd(x)pcs1 = v.compute() # Run on transposition to short-fat array (USV becomes VtSUt)vt, s, ut = da.linalg.svd(x.T)pcs2 = vt.T.compute() # Make sure the PCs are the same (modulo signs)all([ np.allclose(pcs1[i], pcs2[i]) or np.allclose(pcs1[i], -pcs2[i]) for i in range(x.shape[1]) ])> True There are several other references to shapes and axes in the current PCA implementation that would keep this from working as-is, even if it is that simple to get the original components using a transposed input. If there was a way for it to be handled transparently in da.linalg.svd, that would be terrific. It wasn't clear to me in a glance at that code why compute_svd is an option on tsqr but not sfqr. Maybe it would be trivial otherwise? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#731 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AACKZTAZ6WQML53ORZYGTYTSC3JHBANCNFSM4QNJJ7UA> .

[eric-czech]

In any event, I'm +1 if anyone wants to pick this up

Which of these makes the most sense to pick up @mrocklin?

1. Have linalg.svd identify short-fat arrays, transpose them, run the normal tsqr code, and then transpose results (assuming there are no other pitfalls with that)
2. Have some equivalent to tsqr that can operate on short-fat arrays. sfqr sounds like this but on closer inspection it looks like a neglected counterpart to tsqr that has fallen well behind. The docstring still calls it a "quick hack" in any case.

The second option sounds better, but are there many contributors that know the tsqr code well still and could advise on a sfqr implementation / revamp?

I would happily take on #1, and possibly #2 if I could get some substantial help from the start.

[TomAugspurger]

If option 1 works then it sounds like less work. I'm not sure that we have anyone who is really an expert on the tsqr code.

…

On Fri, Aug 28, 2020 at 5:30 AM Eric Czech ***@***.***> wrote: In any event, I'm +1 if anyone wants to pick this up Which of these makes the most sense to pick up @mrocklin <https://github.com/mrocklin>? 1. Have linalg.svd identify short-fat arrays, transpose them, run the normal tsqr code, and then transpose results (assuming there are no other pitfalls with that) 2. Have some equivalent to tsqr that can operate on short-fat arrays. sfqr sounds like this but on closer inspection it looks like a neglected counterpart to tsqr that has fallen well behind. The docstring still calls it a "quick hack" in any case. The second option sounds better, but are there many contributors that know the tsqr code well still and could advise on a sfqr implementation / revamp? I would happily take on #1 <#1>, and possibly #2 <#2> if I could get some substantial help from the start. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#731 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAKAOIUF733TKKOZKLHW2Z3SC6BMDANCNFSM4QNJJ7UA> .

[mrocklin]

+1 on option 1. I think that tsqr should satisfy our needs. And if we can have only one really complicated function rather than two then that would be my preference for maintenance reasons. On Fri, Aug 28, 2020 at 4:15 AM Tom Augspurger <notifications@github.com> wrote:

…

If option 1 works then it sounds like less work. I'm not sure that we have anyone who is really an expert on the tsqr code. On Fri, Aug 28, 2020 at 5:30 AM Eric Czech ***@***.***> wrote: > In any event, I'm +1 if anyone wants to pick this up > > Which of these makes the most sense to pick up @mrocklin > <https://github.com/mrocklin>? > > 1. Have linalg.svd identify short-fat arrays, transpose them, run the > normal tsqr code, and then transpose results (assuming there are no other > pitfalls with that) > 2. Have some equivalent to tsqr that can operate on short-fat arrays. > sfqr sounds like this but on closer inspection it looks like a > neglected counterpart to tsqr that has fallen well behind. The docstring > still calls it a "quick hack" in any case. > > The second option sounds better, but are there many contributors that know > the tsqr code well still and could advise on a sfqr implementation / revamp? > > I would happily take on #1 <#1>, and > possibly #2 <#2> if I could get some > substantial help from the start. > > — > You are receiving this because you were mentioned. > Reply to this email directly, view it on GitHub > <#731 (comment)>, or > unsubscribe > < https://github.com/notifications/unsubscribe-auth/AAKAOIUF733TKKOZKLHW2Z3SC6BMDANCNFSM4QNJJ7UA > > . > — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#731 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AACKZTHCTQVJJS6YG6NSXOLSC6GT7ANCNFSM4QNJJ7UA> .

[eric-czech]

Feel free to close this, otherwise I'll leave it open until I get a chance to explore those changes some more.

[TomAugspurger]

I'm happy to leave this open till we have a resolution.

…

On Fri, Aug 28, 2020 at 9:54 AM Eric Czech ***@***.***> wrote: Feel free to close this, otherwise I'll leave it open until I get a chance to explore those changes some more. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#731 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAKAOIWUCPF55Y6RQXZE7N3SC7AKPANCNFSM4QNJJ7UA> .

[jerowe]

@eric-czech I was just coming here to raise this issue!

I have a PR with some different shaped arrays to demonstrate that dask-ml PCA fails on skinny arrays, which would be mostly all of our arrays. ;-)

https://github.com/pystatgen/sgkit/blob/6dcf446583bc52de5e15bae65f4aa63fcc5bb9b7/sgkit/tests/test_decomposition.py#L239