[Paper Review] Beyond Pham's algorithm for joint diagonalization

Pierre Ablin, J.-F. Cardoso+1 more|arXiv (Cornell University)|Nov 28, 2018

Blind Source Separation Techniques8 references7 citations

TL;DR

This paper proposes a novel quasi-Newton method for approximate joint diagonalization of symmetric matrices using a sparse Hessian approximation that matches the true Hessian when matrices are diagonal. The method achieves faster convergence than Pham’s algorithm, outperforming it by an order of magnitude in numerical experiments on synthetic and real MEG data, with a per-iteration cost of O(np²) and quadratic convergence near optimal solutions.

ABSTRACT

The approximate joint diagonalization of a set of matrices consists in finding a basis in which these matrices are as diagonal as possible. This problem naturally appears in several statistical learning tasks such as blind signal separation. We consider the diagonalization criterion studied in a seminal paper by Pham (2001), and propose a new quasi-Newton method for its optimization. Through numerical experiments on simulated and real datasets, we show that the proposed method outper-forms Pham's algorithm. An open source Python package is released.

Motivation & Objective

• Address the need for faster and more efficient optimization in approximate joint diagonalization, a key task in blind source separation and statistical learning.
• Overcome the high computational cost of Newton-type methods, which require solving O(p⁴) systems, by introducing a sparse Hessian approximation.
• Improve upon Pham’s block coordinate descent algorithm, which, while convergent, is slower in practice, especially on noisy or non-exact diagonalizable sets.
• Develop a method that maintains quadratic convergence near optimal solutions while scaling efficiently with problem size.
• Provide a practical, open-source implementation for use in signal processing and machine learning applications.

Proposed method

• Use relative updates on the invertible matrix B via B ← (I + 𝔼)B to exploit the group structure of the optimization manifold.
• Derive the relative gradient G and relative Hessian 𝒫 from the second-order Taylor expansion of the Pham criterion L(B).
• Approximate the Hessian tensor 𝒫 using a sparse, structured form that depends only on diagonal elements of Dⁱ = B Cⁱ Bᵀ, reducing computation to O(np²).
• Construct a Hessian inverse approximation 𝒫⁺ that enables fast search direction computation via a closed-form solution involving Γ_ab = (1/n)Σᵢ Dⁱ_bb / Dⁱ_aa.
• Use backtracking line search to ensure sufficient decrease in the cost function at each step, maintaining global convergence despite the quasi-Newton approximation.
• Leverage the fact that the Hessian approximation becomes exact when the matrices are exactly jointly diagonalizable, enabling quadratic convergence near the solution.

Experimental results

Research questions

• RQ1 Can a quasi-Newton method be designed for joint diagonalization that achieves faster convergence than Pham’s algorithm while maintaining low per-iteration cost?
• RQ2 Does a sparse Hessian approximation that matches the true Hessian in the jointly diagonalizable case preserve quadratic convergence properties?
• RQ3 Can the proposed method scale efficiently to large p while remaining practical, avoiding the O(p⁶) cost of full Hessian inversion?
• RQ4 How does the method perform on real-world data, such as MEG signals, where exact joint diagonalization is not achievable?
• RQ5 What is the impact of noise on convergence speed and final performance compared to Pham’s algorithm?

Key findings

• The proposed quasi-Newton method outperforms Pham’s algorithm by approximately one order of magnitude in convergence speed on both synthetic and real MEG datasets.
• On the perfectly jointly diagonalizable synthetic dataset (σ = 0), both methods exhibit quadratic convergence, confirming theoretical expectations.
• In the presence of noise (σ = 0.1) and on real MEG data, convergence slows to linear rate, as expected, but the proposed method still converges significantly faster than Pham’s algorithm.
• The Hessian approximation is exact when the set is jointly diagonalizable, ensuring that the method inherits the quadratic convergence rate of Newton’s method in that regime.
• The method achieves a per-iteration cost of O(np²), matching the natural complexity of the problem and enabling scalability to larger p.
• An open-source Python package is released, facilitating reproducibility and integration into downstream applications in signal processing and machine learning.