This reproducible R Markdown analysis was created with workflowr (version 1.7.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(1) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version f421e5c. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were made to the R Markdown (analysis/fastICA_asymmetric.Rmd) and HTML (docs/fastICA_asymmetric.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
Standard fastICA maximizes the log-cosh contrast function, which is equivalent to assuming a symmetric binary (Rademacher) prior on the independent components. For sparse or skewed sources (e.g. indicator variables where the “on” fraction \(p \ll 0.5\)), the expected log-cosh contrast falls below the Gaussian baseline, causing fastICA to actively avoid the true source direction.
This analysis generalizes the framework to an asymmetric Rademacher prior parameterized by its up-probability \(p\). We alternate between:
Updating the projection direction \(w\) via a Newton-like fixed-point step (for fixed \(p\)).
Updating \(p\) by maximizing the marginal log-likelihood (for fixed \(w\)).
At \(p = 0.5\) the algorithm exactly reduces to standard log-cosh fastICA. Two scale anchors are compared:
M-estimator (\(c = s\)): robust but slightly over-inflated variance.
Golden-ratio (\(c \approx 0.618\,s\)): satisfies \(c^2 + cs = s^2\), so the assumed generative variance equals the empirical data variance.
Helpers
prewhiten = function(X, n.comp) {
X = X - rowMeans(X)
sqrt(ncol(X)) * t(svd(X)$v[, 1:n.comp])
}
Implementation
With \(\sigma^2 = cs\), the score function and its derivative are:
Two diagonal Hessian approximations are supported, following the notation in ebproj_newton:
fastICA: \(\tilde H = \overline{M'(x)}\,\mathbf{I}\) — uniform sample weights.
trace: \(\tilde H = c_\text{trace}\,\mathbf{I}\) where \(c_\text{trace} = \tfrac{1}{k}\sum_i M'(x_i)\,S_{ii}\) and \(S_{ii} = \tfrac{1}{n}\|Y_{:i}\|^2\) is the squared distance of sample \(i\) from the origin (with \(\sum_i S_{ii} = k\) for whitened data).
Both reduce to the same Newton fixed-point update structure: \[w \leftarrow \tfrac{1}{n}Y M(x) - \tilde H\,w, \qquad w \leftarrow w / \|w\|\]
fastica_asym_r1 = function(Y, s = 1, anchor = c("M", "golden"),
hess = c("fastICA", "trace"),
tol = 1e-6, max_iter = 500, eps = 0.01,
w_init = NULL) {
anchor = match.arg(anchor)
hess = match.arg(hess)
c = if (anchor == "M") s else s * (sqrt(5)-1)/2
m = nrow(Y); n = ncol(Y)
S_diag = colSums(Y^2) / n # S_ii = ||Y[:,i]||^2 / n (sum = m)
w = if (is.null(w_init)) rnorm(m) else w_init
w = w / sqrt(sum(w^2))
p = 0.5
for (iter in seq_len(max_iter)) {
w_old = w; p_old = p
x = as.vector(t(Y) %*% w)
sc = asym_score(x, p, c, s)
h = if (hess == "trace") sum(sc$Mp * S_diag) / m else mean(sc$Mp)
w = as.vector(Y %*% sc$M) / n - h * w
w = w / sqrt(sum(w^2))
x = as.vector(t(Y) %*% w)
opt = optimize(\(pp) asym_obj_p(pp, x, c, s), c(eps, 1-eps), maximum = TRUE)
p = opt$maximum
if (1 - abs(sum(w * w_old)) < tol && abs(p - p_old) < tol) break
}
list(w = w, p = p, iter = iter, c = c)
}
Standard log-cosh fastICA for comparison:
fastica_r1update = function(X, w) {
w = w / sqrt(sum(w^2))
P = t(X) %*% w
G = tanh(P); G2 = 1 - tanh(P)^2
w = X %*% G - mean(G2) * ncol(X) * w
w / sqrt(sum(w^2))
}
run_seeds_lc = function(Y, S_true, n_seeds = 100, n_iter = 200) {
maxcor = numeric(n_seeds)
for (seed in seq_len(n_seeds)) {
set.seed(seed)
w = rnorm(nrow(Y))
for (i in seq_len(n_iter)) w = fastica_r1update(Y, w)
maxcor[seed] = max(abs(cor(t(S_true), t(Y) %*% w)))
}
maxcor
}
run_seeds_asym = function(Y, S_true, anchor, hess = "fastICA", n_seeds = 100) {
maxcor = numeric(n_seeds); ps = numeric(n_seeds)
for (seed in seq_len(n_seeds)) {
set.seed(seed)
res = fastica_asym_r1(Y, anchor = anchor, hess = hess, w_init = rnorm(nrow(Y)))
maxcor[seed] = max(abs(cor(t(S_true), t(Y) %*% res$w)))
ps[seed] = min(res$p, 1 - res$p) # probability of the rare state
}
list(maxcor = maxcor, p = ps)
}
# Warm-start variant: run log-cosh to convergence, then hand off to asymmetric
run_seeds_asym_warm = function(Y, S_true, anchor, hess = "fastICA", n_seeds = 100,
n_iter_lc = 200) {
maxcor = numeric(n_seeds); ps = numeric(n_seeds)
for (seed in seq_len(n_seeds)) {
set.seed(seed)
w = rnorm(nrow(Y))
for (i in seq_len(n_iter_lc)) w = fastica_r1update(Y, w)
res = fastica_asym_r1(Y, anchor = anchor, hess = hess, w_init = w)
maxcor[seed] = max(abs(cor(t(S_true), t(Y) %*% res$w)))
ps[seed] = min(res$p, 1 - res$p)
}
list(maxcor = maxcor, p = ps)
}
Sanity check: \(p = 0.5\) recovers log-cosh
At \(p = 0.5\), \(s = 1\): \(\Delta(x) = 2x\), \(\pi(x) = (1 + \tanh x)/2\), and \(M(x) = \tanh(x)\) — the standard fastICA score.
x_grid = seq(-3, 3, length.out = 300)
sc05 = asym_score(x_grid, p = 0.5, c = 1, s = 1)
plot(x_grid, sc05$M, type = "l", col = "steelblue", lwd = 2,
xlab = "x", ylab = "M(x)",
main = "Score function at p = 0.5 vs tanh(x)")
lines(x_grid, tanh(x_grid), col = "tomato", lty = 2, lwd = 2)
legend("topleft", c("asymmetric M(x), p = 0.5", "tanh(x)"),
col = c("steelblue", "tomato"), lty = c(1, 2), lwd = 2, bty = "n")
Warning! The custom fig.path you set was ignored by workflowr.
Standard fastICA with log-cosh actively avoids sources with $p < $ about \(0.1\) or \(p > 0.9\): the contrast falls below the Gaussian baseline, so the algorithm prefers noise directions over the true source.
Score functions for varying \(p\)
As \(p\) decreases below 0.5 the score shifts and steepens, penalising the positive tail more heavily — appropriate for sources that are rarely “on”.
pvec2 = c(0.05, 0.1, 0.2, 0.3, 0.5)
cols = c("purple", "tomato", "darkorange", "steelblue", "black")
plot(NULL, xlim = c(-3, 3), ylim = c(-2.5, 2.5),
xlab = "x", ylab = "M(x)",
main = "Asymmetric score functions (c = s = 1)")
for (i in seq_along(pvec2))
lines(x_grid, asym_score(x_grid, pvec2[i], 1, 1)$M, col = cols[i], lwd = 2)
legend("topleft", paste0("p = ", pvec2), col = cols, lwd = 2, bty = "n")
abline(h = 0, lty = 3, col = "grey60")
Warning! The custom fig.path you set was ignored by workflowr.
Log-cosh fails completely. The asymmetric golden-ratio anchor with fastICA Hessian succeeds from random starts (95%); the trace Hessian is compared directly alongside it.
Why does the golden-ratio anchor do better?
The two anchors produce different score functions via the \(c\)-dependent bias term in \(\Delta(x)\). The golden-ratio anchor (\(c = 0.618s\)) satisfies \(c^2 + cs = s^2\), so the total assumed generative variance matches the empirical variance \(s^2\). The M-estimator (\(c = s\)) over-inflates the assumed variance to \(2s^2\), shifting the logistic midpoint and softening the asymmetry penalty.
We can visualise this: at \(p = 0.2\), \(s = 1\), the two anchors produce noticeably different score functions:
x_grid2 = seq(-4, 4, length.out = 400)
sc_M = asym_score(x_grid2, p = 0.2, c = 1, s = 1)
sc_gr = asym_score(x_grid2, p = 0.2, c = (sqrt(5)-1)/2, s = 1)
plot(x_grid2, sc_M$M, type = "l", col = "steelblue", lwd = 2,
xlab = "x", ylab = "M(x)",
main = "Score functions at p = 0.2: M-estimator vs golden-ratio")
lines(x_grid2, sc_gr$M, col = "tomato", lwd = 2)
lines(x_grid2, tanh(x_grid2), col = "grey50", lty = 2, lwd = 1.5)
legend("topleft",
c("M-estimator (c = s)", "golden-ratio (c = 0.618s)", "tanh (p = 0.5)"),
col = c("steelblue", "tomato", "grey50"),
lty = c(1, 1, 2), lwd = c(2, 2, 1.5), bty = "n")
Warning! The custom fig.path you set was ignored by workflowr.
asym golden (warm, trace) mean = 0.995 frac > 0.9 = 0.99 mean_p = 0.447
Asymmetry parameter recovery (\(k = 1\))
Using \(k = 1\) whitening (which isolates each source exactly), we verify that the estimated \(\hat p\) tracks the true sparse fraction. Because of sign ambiguity in the ICA direction, we report \(\min(\hat p,\, 1-\hat p)\), i.e. the probability of the rare state.
Both anchors track the true sparse fraction closely across \(p \in [0.05, 0.5]\).
Summary
The asymmetric fastICA algorithm alternates between a Newton-like fixed-point update for \(w\) (identical to standard fastICA at \(p = 0.5\)) and 1D optimization of \(p\). Two Hessian approximations are compared: fastICA (isotropic, \(\bar{M'(x)}\,\mathbf{I}\)) and trace (weighted by \(S_{ii} = \|Y_{:i}\|^2/n\)). Key findings from 100 random seeds each:
Setting
log-cosh
golden random fastICA
golden random trace
golden warm fastICA
golden warm trace
Symmetric (\(p=0.5\), \(k=9\))
0.99
0.77
0.74
0.99
0.99
9 groups (\(p\approx 0.2\), \(k=9\))
0.00
0.95
0.93
0.65
0.66
Log-cosh fails completely for the 9-groups case because sparse sources have $E[] < $ Gaussian baseline, causing the algorithm to prefer noise directions.
The golden-ratio anchor (\(c \approx 0.618s\)) dramatically outperforms the M-estimator anchor. Its assumed generative variance matches the empirical variance (\(c^2 + cs = s^2\)), giving a better-calibrated asymmetry penalty.
Warm-starting has opposite effects depending on source type:
Symmetric sources: a random start lets \(p\) drift from 0.5 before \(w\) has converged, degrading performance. Warm-starting from log-cosh avoids this and fully recovers high success rates.
Asymmetric sources: log-cosh converges to a wrong direction (noise PC), and the asymmetric method then inherits that bad start. A random start performs better because it can reach the sparse source from a neutral position.
Trace vs fastICA Hessian: the two approximations perform nearly identically here. The fastICA (isotropic) Hessian is marginally better on random starts (95% vs 93% for 9-groups; 77% vs 74% for symmetric), while warm-start results are tied. The trace Hessian down-weights samples with large \(S_{ii}\) (far from the origin) in proportion to \(M'(x_i)\), but for whitened data this provides no consistent advantage.
Asymmetry recovery (\(k=1\) whitening): \(\hat p = \min(p, 1-p)\) correctly tracks the true sparse fraction over the range \([0.05, 0.5]\).
sessionInfo()
R version 4.4.2 (2024-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS 26.5.2
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
locale:
[1] C
time zone: America/Chicago
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] vctrs_0.7.2 cli_3.6.5 knitr_1.51 rlang_1.1.7
[5] xfun_0.56 stringi_1.8.7 otel_0.2.0 promises_1.5.0
[9] jsonlite_2.0.0 glue_1.8.0 workflowr_1.7.2 rprojroot_2.1.1
[13] git2r_0.36.2 htmltools_0.5.9 httpuv_1.6.16 sass_0.4.10
[17] rmarkdown_2.30 tibble_3.3.1 evaluate_1.0.5 jquerylib_0.1.4
[21] fastmap_1.2.0 yaml_2.3.12 lifecycle_1.0.5 whisker_0.4.1
[25] stringr_1.6.0 compiler_4.4.2 fs_1.6.6 pkgconfig_2.0.3
[29] Rcpp_1.1.1 later_1.4.6 digest_0.6.39 R6_2.6.1
[33] pillar_1.11.1 magrittr_2.0.4 bslib_0.10.0 tools_4.4.2
[37] cachem_1.1.0
