PaperlayerView exact revision on arXiv

Source-linked AI summary

Faster independent component analysis by preconditioning with Hessian approximations

Pierre Ablin, Jean-François Cardoso, Alexandre Gramfort

arXiv:1706.08171v3stat.MLstat.AP

TL;DR

Fast and accurate maximum-likelihood ICA is difficult for large real datasets. Picard combines L-BFGS with sparse Hessian-approximation preconditioning to refine curvature information. Extensive experiments show clear runtime gains, especially on real data where the ICA model may not hold exactly.

  • Problem

    Inference time can bottleneck applications using ICA, motivating a fast and accurate method for maximum-likelihood ICA.

  • Method

    Picard preconditions L-BFGS with sparse Hessian approximations, allowing the algorithm to refine those approximations toward the true curvature.

  • Results

    Picard showed clear running-time gains over state-of-the-art methods across synthetic, EEG, fMRI, and image data, converging quickly on real datasets.

  • Takeaways & Limitations

    Picard combines iteration costs similar to simple quasi-Newton methods with substantially better descent directions for the evaluated datasets.

  • Takeaways & Limitations

    The regularization strategy is a trial-and-error heuristic, although it worked uniformly well on the studied datasets.

Abstract

from arXiv · show

Independent Component Analysis (ICA) is a technique for unsupervised exploration of multi-channel data that is widely used in observational sciences. In its classic form, ICA relies on modeling the data as linear mixtures of non-Gaussian independent sources. The maximization of the corresponding likelihood is a challenging problem if it has to be completed quickly and accurately on large sets of real data. We introduce the Preconditioned ICA for Real Data (Picard) algorithm, which is a relative L-BFGS algorithm preconditioned with sparse Hessian approximations. Extensive numerical comparisons to several algorithms of the same class demonstrate the superior performance of the proposed technique, especially on real data, for which the ICA model does not necessarily hold.

1 Introduction

ICA is widely used to estimate statistically independent source signals, but quickly and accurately maximizing its likelihood is difficult, especially when real data depart from the ICA model. The paper addresses this by preconditioning L-BFGS with Hessian approximations.

  • ICA estimates source signals from observations modeled as linear mixtures of statistically independent latent components.
  • Infomax is a maximum-likelihood ICA method, but its stochastic gradient optimization may require hand-tuning, fail to converge, or converge slowly.
  • Curvature-based methods use likelihood second derivatives to accelerate optimization beyond gradient-only approaches.
  • The complete Hessian is costly to evaluate and invert for large datasets, and requires regularization because the ICA likelihood is non-convex.
  • Hessian approximations have low iteration cost but lack sufficient accuracy on real data, where correlated sources can violate the ICA model.
  • Picard learns curvature through L-BFGS while using Hessian approximations as preconditioners.
  • The paper evaluates Picard and related algorithms through numerical experiments on synthetic, EEG, fMRI, and image data.

2 Likelihood and derivatives

The paper formulates maximum-likelihood ICA through an unmixing matrix and derives relative gradients and Hessians for optimization. Sparse, block-diagonal Hessian approximations are cheaper to compute and invert, but require regularization and are used as preconditioners when their accuracy is limited.

  • Likelihood formulation: ICA models observations as X = AS, where A is an unknown mixing matrix and S contains statistically independent source rows.
  • Likelihood formulation: The objective minimizes the negative averaged log-likelihood L(W) with respect to the unmixing matrix W = A^-1.
  • Likelihood formulation: The likelihood optimization is non-convex because permuting the columns of any minimizing W produces another equivalent minimizer.
  • Likelihood formulation: The paper fixes source densities, using the standard Infomax density model in its experiments, while real data are not expected to satisfy the ICA mixture model exactly.
  • Relative derivatives: Relative variations use a first-order relative gradient G and a second-order relative Hessian H in the Taylor expansion of L((I + E)W).
  • Relative derivatives: The relative Hessian has O(N^3) non-zero coefficients, but computing it requires Θ(N^3 × T) operations, motivating cheaper approximations.
  • Hessian approximations: The ˜H2 approximation replaces ˆh_ijl with δ_jlˆh_ij, while ˜H1 further replaces ˆh_ij with ˆh_iˆσ_j^2 for i ≠ j.
  • Hessian approximations: Both approximations are block diagonal, with only two non-zero coefficient positions for each off-diagonal pair and ˜H_ijji = 1.

3 Preconditioned ICA for Real Data

Picard combines L-BFGS curvature learning with sparse Hessian approximations used as a preconditioner. The algorithm retains low-memory recursive updates while incorporating a current curvature estimate into its search direction.

  • Picard algorithm: The preconditioned L-BFGS procedure computes the relative gradient, forms and regularizes a Hessian approximation, then applies the recursive L-BFGS direction calculation.The resulting update is W_k+1 = (I + α_kp_k)W_k.
  • L-BFGS background: L-BFGS stores only the last m Hessian updates and inverts its approximation recursively instead of storing a potentially oversized dense matrix.The memory parameter m controls how many previous updates contribute to the approximation.
  • Picard algorithm: Picard uses Hessian approximations to initialize the recursive L-BFGS formula for computing the search direction.The preconditioner is updated from the current Hessian approximation while L-BFGS retains recent gradient and iterate information.
  • Line search: Picard’s line search selects a step length for each iteration, while Wolfe-based searches require repeated likelihood and gradient evaluations.Backtracking is cheaper per trial because it requires only a likelihood evaluation, but it does not enforce Wolfe conditions.

4 Related work

The related methods use gradients, Hessian information, or Hessian approximations in different ways to optimize the ICA likelihood. Their trade-offs involve stochastic noise, computational cost, approximation accuracy, and regularization.

  • Infomax: Infomax uses stochastic relative-gradient updates on randomly selected mini-batches, accelerating initial progress on large datasets.Mini-batches can begin decreasing the objective after one batch, whereas full-batch methods require a complete data pass before progress starts.
  • Infomax: Plain stochastic gradients require careful learning-rate annealing and may plateau because the full-data gradient does not reach zero.The plateau level is proportional to the step size, which creates a trade-off between progress and stability.
  • Truncated Newton’s method: Truncated Newton approximates H^-1G with conjugate gradients and stops early, avoiding explicit Hessian construction while retaining a useful Newton-like direction.A Hessian-vector product can be computed at roughly gradient-evaluation cost, with complexity Θ(N^2 × T).
  • Truncated Newton’s method: Hessian approximations can precondition conjugate gradients, roughly halving the iterations needed to reach a given error when solving H^-1G.The approach still faces Hessian regularization difficulties because the smallest eigenvalue is not straightforward to obtain without computing H.
  • Quasi-Newton and trust-region methods: Simple quasi-Newton methods replace H with a positive-definite Hessian approximation, while trust-region ICA minimizes a local quadratic model within a trust region.The experiments distinguish simple quasi-Newton H1, simple quasi-Newton H2, and trust region ICA.

5 Experiments

The experiments were designed to compare algorithms fairly under controlled computational conditions. They measured convergence using gradient norms across repeated runs and included several L-BFGS and Hessian-based variants.

  • Experimental protocol: All implementations used the same sample-scaling operations, so convergence-speed differences were attributed to algorithmic behavior rather than differing core routines.The costly likelihood, score, derivative, gradient, Hessian-approximation, and Hessian-free computations were shared.
  • Experimental protocol: Experiments ran on one core of an Intel Core i7-6600U at 2.6 GHz using NumPy with Intel MKL and numexpr.The setup also optimized evaluations of log cosh(y_i(t)/2) and tanh(y_i(t)/2).
  • Evaluation: Convergence was tracked with the gradient infinity norm, and runs stopped after an iteration limit or when the norm reached typically 10^-8.The norm is defined as max_ij|G_ij|.
  • Evaluation: Median convergence curves summarized repeated experiments, with half the runs finishing faster and half slower than the plotted curve.This presentation was chosen to maintain readability across multiple runs.
  • Compared algorithms: The comparison included vanilla L-BFGS, Picard with H1 and H2, simple quasi-Newton, truncated Newton, and trust-region methods.The implementations of the principal algorithms were reimplemented in Python and made available online.

5.3 Simulation study

The simulation study evaluates ICA algorithms on synthetic mixtures spanning ideal independent sources, unrecoverable source types, and mixtures whose identifiability is limited by finite samples.

  • Simulation design: The synthetic setup varied the number of sources, number of samples, and source densities before randomly mixing and preprocessing the signals.Each source contributed independently sampled observations, and experiments were repeated with different random seeds.
  • Experiment A: Experiment A used T = 10000 samples from N = 50 independent super-Gaussian sources.The ICA assumption holds perfectly in this setup, and ψ = tanh(·/2) is appropriate.
  • Experiment B: Experiment B used T = 10000 samples from N = 15 sources spanning Laplace, Gaussian, and sub-Gaussian densities.The Gaussian and sub-Gaussian sources are not recoverable under the stated model and score choice.
  • Experiment C: Experiment C used T = 5000 samples from N = 40 sources with linearly varying mixtures of Gaussian distributions.The mixture is identifiable, but finite samples make the most Gaussian-like sources difficult to separate.

5.4 Experiments on EEG data

The real-data comparisons include EEG and fMRI datasets, evaluated by tracking the infinity norm of the relative gradient over time and iterations. EEG violates the ICA assumptions because of noise and artifacts, while fMRI is benchmarked using a 40 × ≃60000 signal matrix.

  • EEG datasets: The EEG evaluation uses 13 public datasets with n = 71 signals and T ≃75000 samples after down-sampling.The recordings measure electric-potential changes induced by brain activity and contain noise and artifacts.
  • Comparison setup: Figure 3 compares convergence on EEG, fMRI, and image-patch data using the infinity norm of the relative gradient versus time and iterations.Solid lines denote algorithms using approximate Hessians; dashed lines denote standard counterparts.
  • fMRI dataset: The fMRI benchmark uses CanICA-preprocessed ADHD-200 data with N = 40 and T ≃60000.CanICA constructs a signal matrix from several patients before classical ICA is applied.

5.6 Experiments on natural images

The natural-image experiment applies ICA to 30,000 centered and standardized 8 × 8 patches, producing a 64 × 30000 matrix. On these image data, Picard, standard L-BFGS, and truncated Newton converge in seconds, whereas other methods converge slowly.

  • Representation: ICA can treat columns of the mixing matrix W^-1 as dictionary atoms learned from random image patches.The atoms are interpreted as features of the image data.
  • Dataset: The experiment uses 100 natural open-country images and T = 30000 patches of side 8 × 8 pixels, yielding a 64 × 30000 data matrix.Patches are centered and scaled to mean 0 and variance 1 before whitening.
  • Synthetic comparison: On the synthetic experiment where the ICA mixture model holds, second-order algorithms converge in a handful of iterations, with simple quasi-Newton methods fastest.Picard does not significantly improve over the Hessian approximations in this well-specified setting.
  • Synthetic comparison: On simulations with Gaussian signals, first-order methods perform poorly, while Picard and truncated Newton retain quick convergence as Hessian-only methods slow down.The ICA model is not identifiable in these two simulations because of the Gaussian signals.
  • Results: Picard, standard L-BFGS, and truncated Newton converge in a few seconds on images, while the other algorithms show very slow linear convergence.This result contrasts sharply with the behavior of methods relying only on Hessian approximations on the image dataset.
  • Synthetic comparison: Truncated Newton converges in fewer iterations than Picard, but its added cost means the reduced iteration count does not compensate in CPU time.Across the experiments, Picard is faster when comparing elapsed time.

5.8 Complexity comparison of truncated Newton and preconditioned L-BFGS

The complexity comparison models large-real-data workloads where T is much larger than N^2 and focuses on operations scaling with T. Picard requires about 2×tG per descent direction, whereas truncated Newton requires (2 + Ncg)×tG and is therefore not faster in the reported setting.

  • Cost model: The costly operations scale as Θ(N^2 × T), including gradient computation, H̃2 computation, and truncated-Newton Hessian-free products.Operations scaling as Θ(N × T) are assumed negligible for the simplified analysis.
  • Assumptions: The analysis assumes T is large relative to N^2, as in most real-data applications, and neglects operations that do not scale with T.On EEG, operations not scaling with T account for less than 1% of Picard's total timing.
  • Picard cost: Picard finds each descent direction in about 2×tG, matching the simple quasi-Newton method under the stated cost model.It uses the current gradient and Hessian approximations; remaining operations do not scale with T.
  • Truncated Newton cost: Truncated Newton requires (2 + Ncg)×tG per descent direction because each conjugate-gradient inner loop adds a Hessian-free product.Its extra curvature computation increases per-iteration cost relative to Picard.
  • Trade-off: Although truncated Newton converges in about half as many iterations as Picard, it typically needs more than 10 conjugate-gradient iterations for a satisfying Newton direction.Restricting it to Ncg = 2 produces a direction far from Newton's direction and drastically increases iterations.
  • Conclusion: The analysis concludes that truncated Newton, as described, cannot be faster than Picard.Its lower iteration count does not offset the added per-iteration cost in the examined setting.

5.9 Study of the control parameters of Picard

Picard exposes four control parameters: Hessian approximation choice, L-BFGS memory, line-search trials, and regularization. Across the tested synthetic and real signals, m = 7, nls = 10, and λmin = 10^-2 produced uniformly good performance.

  • Parameters: Picard has four control parameters: the Hessian approximation, L-BFGS memory m, line-search tries nls, and regularization constant λmin.The Hessian choice is binary.
  • Hessian approximation: H̃2 is generally a better preconditioner than H̃1, although their performance difference can be small.This preference is based on the reported experiments.
  • Memory: L-BFGS memory has barely any effect for 3 ≤ m ≤ 15; smaller values lack information, while larger values bias the approximation toward older landscape regions.The reported effect concerns the quality of the Hessian approximation.
  • Line search and regularization: Line-search trials affect convergence speed, and the optimal regularization depends on problem difficulty.Despite this dependence, one parameter setting performed uniformly well across all tested signal types.
  • Recommended setting: The experiments used m = 7, nls = 10, and λmin = 10^-2 for synthetic, EEG, fMRI, and image data.These shared settings yielded uniformly good performance across the tested variety of signals.

6 Conclusion

The work targets fast, accurate maximum-likelihood ICA by combining inexpensive Hessian approximations with L-BFGS preconditioning. Picard retains low iteration cost while improving descent directions and showed faster running time across several data types, but its fixed score functions and heuristic regularization remain limitations.

  • Quasi-Newton Hessian approximations cost only twice as much per iteration as gradient descent but can be inaccurate on real datasets.This limits practical convergence speed relative to what second-order methods might suggest.
  • Picard combines cheap Hessian approximations with L-BFGS preconditioning to improve optimization for maximum-likelihood ICA.The approach refines the approximations to better capture the true curvature.
  • Picard has iteration costs similar to simple quasi-Newton methods while providing far better descent directions.
  • Extensive experiments on synthetic, EEG, fMRI, and image data showed clear running-time gains over the state of the art.The comparisons used careful implementations of several literature methods.
  • The algorithm uses fixed score functions, limiting its current scope for recovering a broader class of sources.Extending it to adaptive scores is identified as future work, while preliminary alternating density-estimation steps impaired convergence speed.
  • Picard’s regularization relies on a trial-and-error heuristic, despite working uniformly well on each studied dataset.More informed strategies could account for Hessian eigenvalues being driven by signal statistics.
Loading 1706.08171v3…