Transcript Document

Matthew L. Wright
Institute for Mathematics and its Applications
University of Minnesota
Preliminaries
Let 𝐾 be a simplicial complex, and let 𝔽2 be the two-element field.
𝐶𝑝 (𝐾) denotes the vector space generated by the 𝑝-dimensional simplices of 𝐾.
𝐶𝑝 (𝐾) consists of all 𝒑-chains, which are formal sums 𝑐 =
𝔽2 and 𝜎𝑗 is a 𝑝-simplex in 𝐾.
𝑗 𝑎𝑗 𝜎_𝑗,
where 𝑎𝑗 ∈
The boundary 𝜕(𝜎𝑗 ) is the formal sum of the (𝑝 − 1)-dimensional faces of 𝜎𝑗 .
A boundary has no boundary: 𝜕 ∘ 𝜕 = 𝜕 2 = 0.
The boundary operator connects the vector spaces into a chain complex 𝐶∗ :
⋯
𝐶𝑝+1
𝜕𝑝+1
𝐶𝑝
𝜕𝑝
𝐶𝑝−1
⋯
Preliminaries
The 𝑝-chains that have no boundary are called 𝒑-cycles, and they form a
subspace 𝑍𝑝 of 𝐶𝑝 .
The 𝑝-chains that are the boundary of (𝑝 + 1)-chains are called 𝒑-boundaries,
and they form a subspace 𝐵𝑝 of 𝐶𝑝 .
Since a boundary has no boundary, 𝐵𝑝 ⊆ 𝑍𝑝 .
The quotient 𝐻𝑝 = 𝑍𝑝 /𝐵𝑝 is the 𝒑th homology group of 𝐾.
The dimension of 𝐻𝑝 is the 𝒑th Betti number of 𝐾 and is denoted 𝛽𝑝 .
Computing Betti Numbers
𝐷𝑝 𝑖, 𝑗 = 1 if the 𝑖 th (𝑝 − 1)-simplex is
a face of the 𝑗th 𝑝-simplex
𝐷0 𝑖, 𝑗 = 0 otherwise
Observe:
rank 𝐵𝑝−1 = rank 𝐷𝑝
rank 𝑍𝑝 = 𝑛𝑝 − rank(𝐷𝑝 )
Thus, 𝛽𝑝 = 𝑛𝑝 − rank(𝐷𝑝 ) − rank 𝐷𝑝+1 .
rank(𝑍𝑝 )
𝐷𝑝
rank(𝐶𝑝−1 )
The boundary matrix 𝐷𝑝 records the face
relationship for 𝑝-simplices and (𝑝 − 1)simplices.
rank(𝐷𝑝 )
rank(𝐵𝑝−1 )
Number the 𝑝-simplices of K
consecutively from 1 to 𝑛𝑝 = rank(𝐶𝑝 ).
rank(𝐶𝑝 )
Reduce 𝐷𝑝 to Smith normal form
to obtain the desired ranks.
Computing Persistent Homology
We start with a filtered simplicial complex:
∅ = 𝐾0 ⊂ 𝐾1 ⊂ ⋯ ⊂ 𝐾𝑚 = 𝐾
Step 1: Sort the simplices to get a total ordering compatible with the
filtration.
Step 2: Obtain a boundary matrix 𝐷 with respect to the total order on
simplices.
Step 3: Reduce the matrix using column additions, always respecting the
total order on simplices.
Step 4: Read the persistence pairs to get the barcode.
Example:
∅ = 𝐾0
⊂
𝐾1
⊂
𝐾2
⊂
𝐾3
⊂
𝐾4 = 𝐾
2
Example:
2
4
1
∅ = 𝐾0
⊂
𝐾1
⊂
6
4
1
1
3
𝐾2
⊂
2
4
1
3
5
𝐾3
⊂
7
6
3
5
𝐾4 = 𝐾
Step 1: Sort the simplices to get a total ordering compatible with the filtration
Compatible means:
1. Simplices in each 𝐾𝑖 precede simplices in 𝐾 − 𝐾𝑖
2. The faces of each simplex precede the simplex itself
2
Example:
2
4
1
∅ = 𝐾0
⊂
⊂
1
3
𝐾2
Step 2: Obtain a boundary matrix 𝐷
with respect to the total order
𝐷 𝑖, 𝑗 = 1 if simplex 𝑖 is a
codimension-1 face of simplex 𝑗
𝐷 𝐼, 𝑗 = 0 otherwise
6
4
1
𝐾1
2
⊂
4
5
𝐾3
1
2
3
4
5
6
7
2
3
5
𝐾4 = 𝐾
⊂
1
7
1
3
6
3
4
5
6
7
1 1
1
1
1 1
1
1
1
=𝐷
2
Example:
2
4
1
∅ = 𝐾0
⊂
⊂
6
4
1
𝐾1
2
1
3
𝐾2
⊂
Step 2: Obtain a boundary matrix 𝐷
with respect to the total order
Define 𝑙𝑜𝑤(𝑗) to be the row number
of the lowest non-zero entry in
column 𝑗. If column 𝑗 is empty, then
𝑙𝑜𝑤 𝑗 = 0.
Example: 𝑙𝑜𝑤 2 = 0 and 𝑙𝑜𝑤 5 = 3
4
5
𝐾3
1
2
3
4
5
6
7
2
3
5
𝐾4 = 𝐾
⊂
1
7
1
3
6
3
4
5
6
7
1 1
1
1
1 1
1
1
1
=𝐷
2
Example:
2
4
1
∅ = 𝐾0
⊂
⊂
6
4
1
𝐾1
2
1
3
𝐾2
⊂
Step 3: Reduce the matrix
5
𝐾3
1
for j = 1 to n:
while ∃ j' < j with low(j') = low(j) ≠ 0:
add column j' to column j
3
4
5
6
7
2
3
5
𝐾4 = 𝐾
⊂
1
6
7
1
3
2
ALGORITHM
4
3
4
5
6
7
1 1 1
0
1
0
1
1 0
1
1
1
1
2
Example:
2
4
1
∅ = 𝐾0
⊂
⊂
6
4
1
𝐾1
2
1
3
𝐾2
⊂
Step 4: Read the persistence pairs
If 𝑙𝑜𝑤 𝑗 = 𝑖 > 0, then simplex 𝑗 is negative,
paired with simplex 𝑖.
4
5
𝐾3
1
2
3
If 𝑙𝑜𝑤 𝑗 = 0, then simplex 𝑗 is positive; look
to row 𝑗 for pairing.
4
If there is no 𝑘 with 𝑙𝑜𝑤 𝑘 = 𝑗, then simplex
𝑗 generates homology in 𝐾.
6
5
7
2
3
5
𝐾4 = 𝐾
⊂
1
7
1
3
6
3
4
5
6
7
1 1
1
1
1
1
1
2
Example:
2
4
1
∅ = 𝐾0
⊂
⊂
6
4
1
𝐾1
2
1
3
𝐾2
simplex 1: positive, unpaired
simplex 2: positive, paired with simplex 4
simplex 3: positive, paired with simplex 5
simplex 4: negative, paired with simplex 2
simplex 5: negative, paired with simplex 3
simplex 6: positive, paired with simplex 7
simplex 7: negative, paired with simplex 6
⊂
4
5
𝐾3
1
2
3
4
5
6
7
2
3
5
𝐾4 = 𝐾
⊂
1
7
1
3
6
3
4
5
6
7
1 1
1
1
1
1
1
2
Example:
2
4
1
∅ = 𝐾0
⊂
𝐾1
⊂
6
4
1
1
3
𝐾2
simplex 1: positive, unpaired
simplex 2: positive, paired with simplex 4
simplex 3: positive, paired with simplex 5
simplex 4: negative, paired with simplex 2
simplex 5: negative, paired with simplex 3
simplex 6: positive, paired with simplex 7
simplex 7: negative, paired with simplex 6
⊂
2
4
1
3
5
𝐾3
7
⊂
6
3
5
𝐾4 = 𝐾
bar 1, ∞ in 𝐻0 homology
simplices 2 and 4 enter at the same time,
so the pair does not produce a bar
bar 2,3 in 𝐻0 homology
bar 3,4 in 𝐻1 homology
2
Example:
2
4
1
∅ = 𝐾0
⊂
𝐾1
⊂
1
3
𝐾2
⊂
4
1
3
5
𝐾3
7
⊂
3
5
𝐾4 = 𝐾
simplices 2 and 4 enter at the same time,
so the pair does not produce a bar
𝐻1
𝐻0
bar 2,3 in 𝐻0 homology
2
6
bar 1, ∞ in 𝐻0 homology
Now we can draw the barcode:
1
6
4
1
2
3
4
bar 3,4 in 𝐻1 homology
Data Structures
1
2
3
4
5
6
7
We store the matrix in a columnsparse format.
Store the matrix as an array of
columns, and each column as a linked
list.
The linked list for column 𝑗 stores the
row indexes of the nonzero entries in
column 𝑗, in reverse order so that the
“lowest” index is accessible at the
beginning of the list.
2
3
3
6
1
1
2
5
4
Data Structures
Suppose we have reduced columns
1 through 𝑗 − 1. If 𝑙𝑜𝑤 𝑗 = 𝑖 > 0,
we want to find 𝑗’ < 𝑗 such that
𝑙𝑜𝑤 𝑗 ′ = 𝑙𝑜𝑤(𝑗).
To avoid looking through all
columns 1 through 𝑗 − 1, we use a
low array to store the index of
each nonempty reduced column
by low number.
low
array
1
2
3
4
5
6
7
1
2 4
2
3
3
6
1
1
2
5
3 5
4
5
6
7
4
Data Structures
Suppose we have reduced columns
1 through 𝑗 − 1. If 𝑙𝑜𝑤 𝑗 = 𝑖 > 0,
we want to find 𝑗’ < 𝑗 such that
𝑙𝑜𝑤 𝑗 ′ = 𝑙𝑜𝑤(𝑗).
To avoid looking through all
columns 1 through 𝑗 − 1, we use a
low array to store the index of
each nonempty reduced column
by low number.
low
array
1
2
3
4
5
6
7
1
2 4
2
3
2
6
1
1
1
5
3 5
4
5
6
7
4
Data Structures
Suppose we have reduced columns
1 through 𝑗 − 1. If 𝑙𝑜𝑤 𝑗 = 𝑖 > 0,
we want to find 𝑗’ < 𝑗 such that
𝑙𝑜𝑤 𝑗 ′ = 𝑙𝑜𝑤(𝑗).
To avoid looking through all
columns 1 through 𝑗 − 1, we use a
low array to store the index of
each nonempty reduced column
by low number.
low
array
1
2
3
4
5
6
7
1
2 4
2
3
6
1
1
5
3 5
4
5
6 7
7
4
“Twist” Optimization
Observations:
reduced matrix
If column 𝑗 is reduced and 𝑙𝑜𝑤 𝑗 = 𝑖 > 0, then
simplex 𝑗 is negative, paired with positive simplex 𝑖.
1
Then reducing column 𝑖 will result in a column of zeros.
2
If 𝑑 is the dimension of simplex 𝑗, then the dimension
of simplex 𝑖 is 𝑑 − 1.
Optimization:
Reduce columns corresponding to higher-dimensional
simplices first.
If column 𝑗 is reduced and 𝑙𝑜𝑤 𝑗 = 𝑖 > 0, then set
column 𝑖 to zero.
1
3
4
5
6
7
2
3
4
5
6
7
1 1
1
1
1
1
1
“Twist” Optimization
Example: same boundary matrix as before
unreduced matrix
1
First reduce columns corresponding to 2-simplices.
Column 7 is reduced, and 𝑙𝑜𝑤 7 = 6.
Thus, negative simplex 7 is paired with positive
simplex 6, so column 6 must be zero.
After reducing all columns corresponding to 2simplices, then reduce columns corresponding to
1-simplices.
Now we are done.
1
2
3
4
5
6
7
2
3
4
5
6
7
1 1
1
1
1 1
1
1
1
Runtime Complexity
Let 𝑛 be the total number of simplices, so 𝐷 is an 𝑛 ×
𝑛 matrix.
We must reduce 𝑛 columns.
Reducing any particular column requires at most 𝑛
column additions.
number of column
additions is 𝑂 𝑛2
Each column addition requires at most 𝑛 operations.
Therefore, the runtime is 𝑂 𝑛3 .
However, in practice, the runtime is often better than cubic.
Note that this complexity is just for the matrix reduction, and doesn’t include
building the filtration or boundary matrix.
Generating Cycles
To find the cycles that represent homology classes, we must keep track of
the column additions.
Consider the operation of adding column 𝑗’ to column 𝑗:
• This is the same as adding the chain represented by column 𝑗’ to that
represented by column 𝑗.
• This corresponds to multiplying 𝐷 on
the right by the elementary matrix
that has 1s on the diagonal and in
entry 𝑗 ′ , 𝑗 .
𝑗
1
𝐷
1
1
1
1
𝑗′
1
1
Generating Cycles
Reducing 𝐷 via column operations is the same as multiplying 𝐷 on the
right by matrix 𝑉.
Matrix 𝑉 is a product of elementary matrices and is upper triangular with
1s on the diagonal.
𝑅
=
𝐷
𝑉
Generating Cycles
The columns of 𝑉 indicate the cycles and chains that kill cycles.
If column 𝑗 in 𝑅 is empty, then column 𝑗 in 𝑉 gives the corresponding cycle.
Otherwise, column 𝑗 in 𝑉 gives the chain that kills the cycle corresponding
to 𝑙𝑜𝑤(𝑗).
=
𝑅
𝐷
𝑉
Bit Tree
Replacing the linked lists (for columns) by a more clever data structure can
make the algorithm faster, though the complexity is still 𝑂 𝑛3 .
Bauer et. al. propose a bit tree: a 64-ary tree, stored as an array.
1
…
…
2
3
…
64
…
…
…
…
…
…
…
Each block has 64 bits. Bit 𝑖 indicates whether the 𝑖 th subtree is non-empty.
The leaves of the tree encode the nonzero entries of a column.
Bit Tree
1
…
…
2
3
…
64
…
…
…
…
…
…
…
The bit tree supports nearly constant time insertion, deletion, and lookup.
The “active” column is converted to a bit tree, then columns are added to
it, and then it is converted back to a sparse structure.
Bauer et. al. recorded significant speed improvements when using the bit
tree for column additions.
Available Software
JavaPlex: http://git.appliedtopology.org/javaplex/
Dionysus: http://www.mrzv.org/software/dionysus/
Perseus: http://www.sas.upenn.edu/~vnanda/perseus/
Persistent Homology Algorithms Toolbox (PHAT):
https://code.google.com/p/phat/
References
Ulrich Bauer, et. al. “PHAT – Persistent Homology Algorithms Toolbox.”
http://phat.googlecode.com
Herbert Edelsbrunner and John Harer. “Persistent homology: a survey.” in
Surveys on discrete and computational geometry: twenty years later. AMS
(2008).
Afra Zomorodian and Gunnar Carlsson. “Computing persistent homology.”
Discrete and Computational Geometry. Vol. 33, no. 2 (2005), p. 249 – 274.