MAGIC Algorithm: Mathematical Foundations

Step 1: Dimensionality Reduction (Optional)

For computational efficiency, MAGIC first reduces dimensionality using PCA:

\[X_{pca} = X \cdot V_k\]

where \(V_k\) contains the top \(k\) principal components (default: 100).

Step 2: k-Nearest Neighbor Graph

For each cell \(i\), identify its \(k\) nearest neighbors based on Euclidean distance in PCA space:

\[\mathcal{N}_k(i) = \{j : d(x_i, x_j) \leq d_k(x_i)\}\]

where \(d_k(x_i)\) is the distance to the \(k\)-th nearest neighbor.

Step 3: α-Decaying Kernel

The affinity between cells is computed using an α-decaying kernel with adaptive bandwidth:

\[K(x_i, x_j) = \exp\left(-\left(\frac{\|x_i - x_j\|_2}{\varepsilon_i}\right)^\alpha\right)\]

where:

• \(\varepsilon_i = d_k(x_i)\) is the adaptive bandwidth (distance to \(k\)-th neighbor)
• \(\alpha\) controls kernel sharpness (default: 1)

Why adaptive bandwidth?

• Accounts for varying local density
• Cells in dense regions have smaller bandwidth
• Cells in sparse regions have larger bandwidth
• Preserves rare cell populations

# Visualize α-decaying kernel
x <- seq(0, 3, length.out = 100)

par(mfrow = c(1, 2))

# Effect of distance
plot(x, exp(-x), type = "l", lwd = 2, col = "#e74c3c",
     xlab = "Normalized Distance (d/ε)", ylab = "Affinity K(x,y)",
     main = "α-Decaying Kernel (α = 1)")
abline(h = 0.5, lty = 2, col = "gray")
abline(v = log(2), lty = 2, col = "gray")

# Effect of α
plot(x, exp(-x^0.5), type = "l", lwd = 2, col = "#3498db",
     xlab = "Normalized Distance (d/ε)", ylab = "Affinity K(x,y)",
     main = "Effect of α Parameter", ylim = c(0, 1))
lines(x, exp(-x^1), lwd = 2, col = "#e74c3c")
lines(x, exp(-x^2), lwd = 2, col = "#2ecc71")
legend("topright", legend = c("α = 0.5", "α = 1", "α = 2"),
       col = c("#3498db", "#e74c3c", "#2ecc71"), lwd = 2)

Step 4: Markov Normalization

The kernel matrix is row-normalized to create a Markov transition matrix:

\[P = D^{-1}K\]

where \(D\) is the diagonal degree matrix with \(D_{ii} = \sum_j K_{ij}\).

Properties of P:

• Row-stochastic: \(\sum_j P_{ij} = 1\)
• Represents transition probabilities in a random walk
• Eigenvalues in \([-1, 1]\)

Step 5: Diffusion (Powering)

The diffusion operator is powered \(t\) times:

\[X_{imputed} = P^t X\]

Each power of \(P\) propagates information one step along the graph:

• \(t = 1\): Local averaging with immediate neighbors
• \(t = 2\): Information from 2-hop neighbors
• \(t = t\): Information from \(t\)-hop neighborhood

# Simulate diffusion effect
set.seed(42)
n <- 100
x <- c(rnorm(50, -2, 0.5), rnorm(50, 2, 0.5))
y <- c(rnorm(50, 0, 0.5), rnorm(50, 0, 0.5))

# Add noise
x_noisy <- x + rnorm(n, 0, 0.3)
y_noisy <- y + rnorm(n, 0, 0.3)

par(mfrow = c(1, 2))
plot(x_noisy, y_noisy, pch = 16, col = adjustcolor("#3498db", 0.6),
     main = "Before Diffusion (Noisy)", xlab = "Gene 1", ylab = "Gene 2")

# Simple smoothing simulation
x_smooth <- 0.7 * x + 0.3 * x_noisy
y_smooth <- 0.7 * y + 0.3 * y_noisy
plot(x_smooth, y_smooth, pch = 16, col = adjustcolor("#e74c3c", 0.6),
     main = "After Diffusion (Denoised)", xlab = "Gene 1", ylab = "Gene 2")

Procrustes Disparity

Given two matrices \(X\) and \(Y\), Procrustes analysis finds the optimal rotation \(R\) that minimizes:

\[d_{Procrustes}(X, Y) = \|X - YR\|_F^2\]

after centering and normalizing both matrices.

Convergence Criterion

MAGIC iteratively applies diffusion and monitors:

\[\Delta_t = d_{Procrustes}(X_t, X_{t+1})\]

The algorithm stops when \(\Delta_t < \theta\) (default: 0.001), indicating convergence.

library(MAGICR)
data(magic_testdata)

# Run with automatic t selection
result <- magic(magic_testdata, t = "auto", t_max = 15, verbose = TRUE)

cat("\nOptimal t selected:", result$params$t, "\n")
#> 
#> Optimal t selected: 15

Exact Solver

Computes imputation in full gene space:

\[X_{imputed} = P^t X\]

Advantages: - Preserves all gene-gene relationships - No approximation error

Disadvantages: - Memory intensive for large gene sets - Slower computation

Approximate Solver

Operates in PCA space and projects back:

\[X_{imputed} = P^t X_{pca} \cdot V_k^T\]

Advantages: - Much faster for large datasets - Lower memory footprint

Disadvantages: - Some information loss from PCA projection

Parameter Selection Guidelines

When to Use MAGIC

✅ Good use cases: - Visualizing gene-gene relationships - Trajectory analysis preprocessing - Recovering dropout in sparse data

❌ Avoid when: - Performing differential expression (use raw counts) - Data is already dense - Preserving exact count values is critical