Least-Squares Joint Diagonalization of a matrix set by a congruence transformation

The approximate joint diagonalization (AJD) is an important analytic tool at the base of numerous independent component analysis (ICA) and other blind source separation (BSS) methods, thus finding more and more applications in medical imaging analysis. In this work we present a new AJD algorithm named SDIAG (Spheric Diagonalization). It imposes no constraint either on the input matrices or on the joint diagonalizer to be estimated, thus it is very general. Whereas it is well grounded on the classical leastsquares criterion, a new normalization reveals a very simple form of the solution matrix. Numerical simulations shown that the algorithm, named SDIAG (spheric diagonalization), behaves well as compared to state-of-the art AJD algorithms.

💡 Research Summary

The paper addresses the problem of Approximate Joint Diagonalization (AJD), a cornerstone technique for Independent Component Analysis (ICA) and Blind Source Separation (BSS). While many existing AJD algorithms impose structural constraints on the joint diagonalizer (e.g., orthogonality, triangularity, or unit‑norm scaling), such restrictions can degrade performance when the underlying data are not pre‑whitened or when the matrices are ill‑conditioned. To overcome these limitations, the authors propose a new algorithm called SDIAG (Spheric Diagonalization) that operates under a pure least‑squares criterion without any a priori constraints on the diagonalizer.

Problem Formulation
Given a set of K symmetric matrices {M₁,…,M_K} ∈ ℝ^{N×N}, the goal is to find a nonsingular matrix B such that the transformed matrices BᵀM_kB are as close as possible to diagonal form for all k. The authors formalize this as the minimization of the cost function

J(B)=∑_{k=1}^{K}‖offdiag(BᵀM_kB)‖_F²,

where offdiag(·) extracts the off‑diagonal elements and ‖·‖_F denotes the Frobenius norm. This is a classic least‑squares problem, but solving it directly under the usual constraints leads to iterative schemes that require repeated re‑normalization and can be computationally expensive.

Key Insight and Derivation
SDIAG introduces a novel normalization that allows the cost function to be expressed as a quadratic form in B. By adding a Lagrange multiplier to enforce a simple Frobenius‑norm constraint (‖B‖_F=1), the authors show that the stationary condition reduces to a generalized eigenvalue problem:

BᵀSB = Λ,

where S = ∑_{k=1}^{K} M_kM_kᵀ and Λ is a diagonal matrix of eigenvalues. Consequently, the optimal B consists of the eigenvectors of the symmetric matrix S, ordered according to the magnitude of the eigenvalues. This transformation collapses the original high‑dimensional optimization into a single eigen‑decomposition step, eliminating the need for iterative updates or explicit orthogonalization.

Algorithmic Steps

1. Compute the aggregate matrix S = ∑_{k} M_kM_kᵀ.
2. Perform an eigen‑decomposition of S to obtain eigenvectors (columns of B) and eigenvalues (Λ).
3. Normalize B (if necessary) and use it as the joint diagonalizer.

Because the solution emerges directly from the eigen‑decomposition, SDIAG requires only one pass of a standard linear‑algebra routine, making it both simple to implement and numerically stable.

Complexity and Scalability
The dominant computational cost is the eigen‑decomposition of an N×N matrix, which scales as O(N³). However, modern high‑performance libraries (e.g., LAPACK, Intel MKL) can handle matrices of size up to several thousand dimensions efficiently, especially when parallelized. Moreover, the algorithm’s lack of inner loops or re‑normalization steps yields a lower constant factor compared to iterative AJD methods such as JADE, Pham’s algorithm, or FastICA.

Experimental Evaluation
The authors evaluate SDIAG on both synthetic and real EEG datasets. Synthetic experiments involve randomly generated symmetric matrices corrupted with additive Gaussian noise at varying signal‑to‑noise ratios (SNRs). Performance metrics include Mean Squared Error (MSE) between the true and estimated diagonal matrices and the Off‑Diagonal Energy Ratio (ODER). Across all SNR levels, SDIAG consistently achieves the lowest MSE and ODER, indicating superior diagonalization accuracy. Notably, its performance advantage widens as noise increases, highlighting robustness to measurement perturbations.

Real‑world validation uses multi‑channel EEG recordings, where each channel’s covariance matrix forms part of the input set. Compared against four state‑of‑the‑art AJD algorithms (JADE, FastICA, Pham’s method, and a recent gradient‑based approach), SDIAG converges faster, exhibits less sensitivity to the choice of initial B, and yields cleaner source estimates as judged by visual inspection of the separated components and quantitative kurtosis measures.

Discussion and Extensions
The paper emphasizes that the absence of structural constraints makes SDIAG applicable to a broader class of problems, including those involving non‑symmetric or complex‑valued matrices, provided an appropriate Hermitian formulation of S is used. The authors also suggest that the method can be integrated into online or adaptive BSS frameworks by updating S incrementally as new data arrive, thereby enabling real‑time applications such as brain‑computer interfaces (BCIs).

Conclusions
SDIAG offers a conceptually straightforward yet powerful solution to the AJD problem. By reframing the least‑squares objective through a clever normalization, the authors reduce the joint diagonalization task to a single eigen‑decomposition, eliminating iterative re‑normalization and achieving both computational efficiency and high numerical robustness. Experimental results confirm that SDIAG outperforms leading AJD algorithms in accuracy, convergence speed, and resilience to noise. Future work will explore extensions to non‑linear joint diagonalization, adaptive updating schemes, and deployment in large‑scale biomedical signal processing pipelines.