bd.approx: Two Block Diagonal Matrix Approximation.

Description

Finds the best two block diagonal matrix approximation of a similarity matrix according to some distance (linkage) function as described in Bauer (202X). By default, candidate splits are determined by sparse approximations of the first eigenvectors of the similarity matrix. Alternatively, all unique two-way partitions can be evaluated by exhaustive search.

Usage

bd.approx(
  S,
  linkage = "average",
  q = 1,
  h.power = 2,
  balance,
  max.iter,
  search = "hcsvd",
  exhaustive.max.p = 20
)

Value

A list containing:

B: The best two block diagonal matrix approximation.
BD: The best two block diagonal matrix approximation permuted to a block diagonal shape: \(BD = P B t(P)\).
P: The permutation matrix \(P\): \(BD = P B t(P)\).
clustering: The clustering vector as an integer vector of length \(p\), which gives for each component the number 1 or 2 of the cluster/split to which it belongs.
split: A list containing the two splits.
distance: The approximation error (distance) according to the selected linkage function.
search: The search strategy used for the split.
n.candidates: The number of candidate evaluations performed.

Arguments

S: A scaled \(p\)x\(p\) similarity matrix. For example, this may be a correlation matrix.
linkage: The linkage function to be used. This should be one of "average", "ncut", "rcut", "RV" (for RV-coefficient), or "single".
q: Number of sparse eigenvectors to be used whenever the HC-SVD search is applied. This should be either a numeric value between zero and one to indicate percentages, or "Kaiser". For a numerical value between zero and one, the number of sparse eigenvectors is determined as the corresponding share of the total number of eigenvectors, capped at \(p-1\). E.g., q = 1 uses all \(p-1\) available sparse eigenvectors. The Kaiser rule uses the number of eigenvalues of the Hadamard-transformed similarity matrix that are at least one, plus two additional eigenvectors, capped at \(p-1\).
h.power: h-th Hadamard power of S used whenever sparse-eigenvector candidate generation is applied. This should be a positive integer. Larger powers shrink similarities with magnitude smaller than one towards zero and can increase the contrast between within- and between-block similarities, although their effect depends on the similarity structure and clustering objective.
balance: Minimum proportion of the smaller block when splitting into two. Must be a numeric value in \((0, 0.5]\). For example, balance = 0.5 enforces an exact 50:50 split, while balance = 0.2 allows splits as unbalanced as 20:80, with more balanced splits such as 30:70 or 40:60 also permitted. If an exact split is not possible (e.g., balance = 0.5 when \(p = 9\)), the closest integer partition is used (e.g., 4 and 5 per block).
max.iter: How many iterations should be performed for computing the sparse eigenvectors when the HC-SVD search is used. Default is 500.
search: Search strategy. "hcsvd" uses sparse-eigenvector candidate generation, "exhaustive" evaluates all unique two-way partitions for the selected linkage function only, and "auto" uses exhaustive search when \(p\) is at most exhaustive.max.p and HC-SVD otherwise.
exhaustive.max.p: Largest problem size for which exhaustive search is used when search = "auto". Default is 20. The exhaustive C++ implementation is restricted to at most 30 objects.

Details

The sparse eigenvectors are computed using the method of Shen and Huang (2008). The method is implemented by Baglama, Reichel, and Lewis in ssvd (irlba). Here, we use a Rcpp/RcppArmadillo implementation based on ssvd with slight modifications to suit our method and for faster performance.

References

Bauer, J.O. (202X). Divisive hierarchical clustering using block diagonal matrix approximations. Working paper.

Shen, H. and Huang, J.Z. (2008). Sparse principal component analysis via regularized low rank matrix approximation, J. Multivar. Anal. 99, 1015–1034.

Examples

Run this code

#We give a trivial example for a block diagonal matrix perturbed by
#noise, for adapting clustering objectives of spectral clustering,
#and for balanced clustering.

# \donttest{
### TOY EXAMPLE

A <- matrix(c(2,1,1,3), 2, 2)  # 2x2 block
B <- matrix(c(5,4,4,6), 2, 2)  # 2x2 block

# Create a 5x5 zero matrix and insert blocks at right positions.
M <- matrix(0, 4, 4)
M[1:2, 1:2] <- A
M[3:4, 3:4] <- B

M.tilde <- M + matrix(rnorm(4^2, 0, 0.2), 4, 4)

#Construct a similaritiy matrix with same block structure
S <- cov2cor(t(M.tilde) %*% M.tilde)
bd <- bd.approx(S)

#Block diagonal approximation:
bd$B

#We can also permute the block diagonal shape:
S2 <- S[c(1, 3, 2, 4), c(1, 3, 2, 4)]
bd2 <- bd.approx(S2)

#bd2$B gives us again the block diagonal approximation
bd2$B

#And bd2$BD gives us the block diagonal approximation permuted to
#block diagonal shape
bd2$BD


### ADAPTING CLUSTERING OBJECTIVES

#We will use the USArrests example (see ?hcsvd).
data("USArrests")
USArrests["Maryland", "UrbanPop"] <- 76.6
D <- as.matrix(dist(USArrests))
S <- 1 - D / max(D)

#We compute k = 2 clusters adapting the objective of spectral clustering
#with the ratio cut.
bd.approx(S, linkage = "rcut")


### BALANCED CLUSTERING

#We can also enforce balanced clustering, such as two clusters of equal
#size (50:50). We will do this for the USArrests example from above.
bd.approx(S, linkage = "rcut", balance = 0.5)
# }

Run the code above in your browser using DataLab