Tiny VoltPaper: A Flag Decomposition for Hierarchical Datasets
This article explains the key idea in the paper A Flag Decomposition for Hierarchical Datasets by Mankovich et al. It allows one to decompose (or factorize) a matrix A into a product of matrices QR which tries to preserve the hierarchical structure of the data.
0. Problem Statement
Hierarchical data can be represented by a list of column vectors , where . We can also represent the extra columns that are added by (assuming a column as an element in the ordered set of columns) with . In the example below, three such matrices are shown.
Note that is the entire data matrix. We want to find an orthogonal basis for each subspace such that the hierarchical structure is preserved.
1. Key idea
We want to create a hierarchy of subspaces that respects the hierarchical structure of the data. Such a subspace can be represented by an orthogonal matrix . Specifically, for each , we want to find a set of orthonormal vectors that preserve the structure of (in a way that will be defined shortly) as much as possible. We also have for , meaning the subspaces are orthogonal to each other.
Once we obtain these matrices, here is how the data can be represented:
where k is the number of levels in the hierarchy. Here is how one can interpret each term in the sum above:
1.1 Finding the matrices
The expression can be seen as finding a projection matrix that best approximates in the subspace spanned by . This is essentially solving a least-squares problem. The solution to it is given by the first columns of the left singular vectors of .
The expression can be seen as finding a projection matrix that best approximates in the subspace spanned by , after removing the component already explained by .
The expression can be seen as finding a projection matrix that best approximates in the subspace spanned by , after removing the component already explained by and .
This pattern continues for all where .
2. Decomposition
Now that we have computed all the matrices, we can write the data matrix as:
The right matrix can be re-written as a block upper triangular matrix where each entry . Note that these arguments hold true for any .
3. Algorithm
The following observation will help us implement the algorithm efficiently:
The implementation below is largely a modification of the official implementation of the algorithm:
def flag_decomposition(data, Aset, flag_type):
n, p = data.shape
k = len(flag_type)
ms = [flag_type[0]] + [flag_type[i]-flag_type[i-1] for i in range(1,k)]
R = np.zeros((sum(ms), p))
Qs = []
B_indices = [np.array(Aset[0])]+[np.setdiff1d(Aset[i],Aset[i-1])for i in range(1,len(Aset))]
Bs = [data[:,Bset_i] for Bset_i in B_indices]
for i in range(k):
for j in range(i,k):
if j == i:
Ui,_,_ = np.linalg.svd(Bs[i], full_matrices=False)
Ui = Ui[:,:ms[i]]
Qs.append(Ui)
R_block = Ui.T @ Bs[i]
R[sum(ms[:i]):sum(ms[:i+1]), B_indices[j]] = R_block
else:
Ui = Qs[i]
R_block = Ui.T @ Bs[j]
R[sum(ms[:i]):sum(ms[:i+1]), B_indices[j]] = R_block
Bs[j] = Bs[j] - Ui @ R_block
return np.concatenate(Qs, axis=1), R