mrzv / dionysus

Library for computing persistent homology
http://mrzv.org/software/dionysus2
Other
144 stars 31 forks source link

Vines question #45

Closed peekxc closed 3 years ago

peekxc commented 3 years ago

This is actually a question regarding the source in Dionysus 1.

Around Line ~120 of dynamic-persistence.hpp looks like:

rLog(rlTranspositions, "l cycle: %s", l->cycle.tostring(outmap).c_str());
        if (!(l->cycle.contains(index(i_prev))))
        {
            // Case 1.2
            rLog(rlTranspositions, "k is in l: %d", (bool) l->trail.contains(index(k)));       // if true, a special update would be needed to maintain lazy decomposition
            swap(i_prev, i);
            rLog(rlTranspositions, "Case 1.2");
            Count(cTranspositionCase12);
            return false;
        }

This is part of the updating the vineyard for Case 1.

I was wondering what the comment "if true, a special update would be needed to maintain lazy decomposition" meant?

mrzv commented 3 years ago

You deserve an award for code archaeology. I don't really remember what this was about — that line was committed in 2009, and it was unrelated to the main part of that commit — but let's see if we can make a guess.

The first question is what did I mean by "lazy decomposition"? An obvious guess would be that it's what ELZ algorithm does by default: reduce the column to the point where you find the pivot, but no further. I believe trails store the rows of matrix U, which means that if the row of l has k in it (the logged condition), then after the transposition the implied column operation would not occur in the original ELZ algorithm/lazy decomposition. So if one wants to maintain that, then an extra update is necessary: zero out U[l,k] by adding R[l] to R[k] (mod 2).

Does this make sense?

peekxc commented 3 years ago

I see, yes I was wondering if lazy meant something like the decomposition is maintained in the sense that R is reduced, U is upper triangular, but you may need additional zero-ing operations to recover the boundary matrix.

I have a situation where I'm trying to maintain an R = DV decomposition instead of an D = RU one. I was following alongside your thesis, but strangely may have found a situation where the operations maintaining D = RU works fine (in the sense that I can permute my original D -> PDP and match it exactly with the one produced by the updated RU), but I'm encountering what may be a corner case (or a mistake on my end), because performing the same (column/column) operations on the R and V matrices isn't always matching the (column/row) operations on the R and U matrices.

More specifically, I have an R = DV decomposition that I believe to be correct, but after I transpose two positive simplices and look at Case 1 my particular matrix has R[i,l]= 0 (so then I permute as usual w/ Case 1.2). I get a decomposition where R is indeed reduced, V is indeed upper-triangular, but when I invert with RV^{-1} I get a matrix where one of the boundary chains is incorrect. So I was looking at your special cases to see if you had an extra condition lying somewhere. Interestingly, though, maintaining the RU directly worked fine

mrzv commented 3 years ago

I see, yes I was wondering if lazy meant something like the decomposition is maintained in the sense that R is reduced, U is upper triangular, but you may need additional zero-ing operations to recover the boundary matrix.

What I called lazy in my comment is unrelated to recovering the boundary matrix. The invariant that D = RU is always maintained, but there are many different decompositions that satisfy it. The point is that it's a bit extra work to maintain precisely the decomposition that ELZ algorithm would compute.

More specifically, I have an R = DV decomposition that I believe to be correct, but after I transpose two positive simplices and look at Case 1 my particular matrix has R[i,l]= 0 (so then I permute as usual w/ Case 1.2). I get a decomposition where R is indeed reduced, V is indeed upper-triangular, but when I invert with RV^{-1} I get a matrix where one of the boundary chains is incorrect. So I was looking at your special cases to see if you had an extra condition lying somewhere.

That doesn't sound related to this condition. It sounds like a bug in the code/algorithm. Are you certain that V is upper-triangular? When maintaining V, you have to worry about V[i,i+1]. Also, does R = DV after the transposition? That seems like the first and more natural condition to check.

mrzv commented 3 years ago

Maybe the problem is with incorrect/misleading phrasing in my thesis? When it describes Case 1, and says that we set U[i,i+1] = V[i,i+1] = 0, for U[i,i+1] you can just set it. For V[i,i+1], you have to add V[i] to V[i+1], otherwise you break the invariant.

peekxc commented 3 years ago

Ahh, adding V[i] to V[i+1] in Case 1 and 4 fixed it! So subtle, I should've realized that before. Thank you @mrzv!

mrzv commented 3 years ago

Glad to hear it worked out. Where are you implementing this? I would very much like to get vineyards back in Dionysus 2. Any interest in creating a pull request?

peekxc commented 3 years ago

I have a package I'm working on for a much larger project which involves vineyards at some level. It's private for now because it's hyper developmental still, but as with all my projects I will eventually make it open source :)

The code is all R/C++ for now, which is sort of the language domain I stick with, so may not be relevant to Dionysus

mrzv commented 3 years ago

Got it. Best of luck with your project.