Principal Component Analysis (PCA) is a foundational technique in machine learning for reducing the dimensionality of high-dimensional datasets. However, PCA can lead to biased representations that disadvantage certain subgroups within the data. To address this issue, a Fair PCA (FPCA) model was introduced to equalize the reconstruction loss between subgroups, but the existing semidefinite relaxation (SDR) based approach is computationally expensive even for a suboptimal solution. Although several alternative FPCA variants have been developed to improve efficiency, they often shift attention away from equalizing the reconstruction loss – the central goal of FPCA. In this paper, we identify a hidden convexity in FPCA and introduce a new algorithm that solves the resulting convex optimization via an eigenvalue optimization. Our approach achieves the desired fairness in reconstruction loss without sacrificing performance. Experiments on real-world datasets show that the proposed FPCA algorithm is approximately 8× faster than the SDR-based algorithm while being at most 85

Zermelo's algorithm is a classical method for computing the maximum likelihood estimator in the Bradley–Terry (BT) model, but its convergence can be slow in practice. To accelerate computation, Newman introduced a family of Zermelo-type fixed-point iterations parameterized by α, with Zermelo's algorithm recovered at α=1. Empirical evidence suggests that the choice α=0 often converges substantially faster, making it a promising alternative, yet the mechanism underlying this acceleration remains elusive. This paper provides theoretical insight into this phenomenon through a systematic local convergence analysis. We derive closed-form expressions for local convergence factors under synchronous and asynchronous updates and analyze their dependence on α via spectral analysis of the associated Jacobian matrices. For synchronous updates, we show that the algorithm may fail to converge when α<1, and its local convergence factor is quasi-convex in α under the population BT model. In contrast, asynchronous updates are always locally convergent, and their local convergence factor is provably monotonically increasing in α under the population BT model of consistently ordered bipartite comparison graphs, establishing the optimality of α=0 in this setting. We further establish asymptotic approximation results for the population convergence factors under the BT model, justifying their practical relevance. Numerical experiments on synthetic and real-world datasets confirm the theory. Our analysis complements existing convergence results and shows that the acceleration of α=0 arises not only from the parameter choice but, more importantly, from the use of asynchronous updates.

Matrix joint block-diagonalization (jbd) frequently arises from diverse applications such as independent component analysis, blind source separation, and common principal component analysis (CPCA), among others. Particularly, CPCA aims at joint diagonalization, i.e., each block size being 1-by-1. This paper is concerned with principal joint block-diagonalization (pjbd), which aims to achieve two goals: 1) partial joint block-diagonalization, and 2) identification of dominant common block-diagonal parts for all involved matrices. This is in contrast to most existing methods, especially the popular ones based on Givens rotation, which focus on full joint diagonalization and quickly become impractical for matrices of even moderate size (300-by-300 or larger). An NPDo approach, directly aiming at dominant common block-diagonal parts, is proposed and it is built on the nonlinear polar decomposition with orthonormal polar factor dependency that characterizes the solutions of the optimization problem designed to achieve pjbd, and it is shown the associated SCF iteration is globally convergent to a stationary point while the objective function increases monotonically during the iterative process. Numerical experiments, including practical applications to multi-view subspace clustering and independent subspace analysis, are presented to illustrate the effectiveness of the NPDo approach and its superiority to Givens rotation-based methods.

We address the problem of computing the eigenvalue backward error of the Rosenbrock system matrix under various types of block perturbations. We establish novel characterizations of these backward errors using a class of minimization problems involving the sum of two generalized Rayleigh quotients (SRQ2). For computational purposes and analysis, we reformulate such optimization problems as minimization of a rational function over the joint numerical range of three Hermitian matrices. This reformulation eliminates certain local minimizers of the original SRQ2 minimization and allows for convenient visualization of the solution. Furthermore, by exploiting the convexity within the joint numerical range, we derive a characterization of the optimal solution using a nonlinear eigenvalue problem with eigenvector dependency (NEPv). The NEPv characterization enables a more efficient solution of the SRQ2 minimization compared to traditional optimization techniques. Our numerical experiments demonstrate the benefits and effectiveness of the NEPv approach for SRQ2 minimization in computing eigenvalue backward errors of Rosenbrock systems.

Principal Component Analysis (PCA) is a foundational technique in machine learning for dimensionality reduction of high-dimensional datasets. However, PCA could lead to biased outcomes that disadvantage certain subgroups of the underlying datasets. To address the bias issue, a Fair PCA (FPCA) model was introduced by Samadi et al. (2018) for equalizing the reconstruction loss between subgroups. The semidefinite relaxation (SDR) based approach proposed by Samadi et al. (2018) is computationally expensive even for suboptimal solutions. To improve efficiency, several alternative variants of the FPCA model have been developed. These variants often shift the focus away from equalizing the reconstruction loss. In this paper, we identify a hidden convexity in the FPCA model and introduce an algorithm for convex optimization via eigenvalue optimization. Our approach achieves the desired fairness in reconstruction loss without sacrificing performance. As demonstrated in real-world datasets, the proposed FPCA algorithm runs 8× faster than the SDR-based algorithm, and only at most 85

We consider a class of eigenvector-dependent nonlinear eigenvalue problems (NEPv) without the unitary invariance property. Those NEPv commonly arise as the first-order optimality conditions of a particular type of optimization problems over the Stiefel manifold, and previously, special cases have been studied in the literature. Two necessary conditions, a definiteness condition and a rank-preserving condition, on an eigenbasis matrix of the NEPv that is a global optimizer of the associated optimization problem are revealed, where the definiteness condition has been known for the special cases previously investigated. We show that, locally close to the eigenbasis matrix satisfying both necessary conditions, the NEPv can be reformulated as a unitarily invariant NEPv, the so-called aligned NEPv, through a basis alignment operation - in other words, the NEPv is locally unitarily invariantizable. Numerically, the NEPv is naturally solved by a self-consistent field (SCF)-type iteration. By exploiting the differentiability of the coefficient matrix of the aligned NEPv, we establish a closed-form local convergence rate for the SCFtype iteration and analyze its level-shifted variant. Numerical experiments confirm our theoretical results.

This paper concerns a class of monotone eigenvalue problems with eigenvector nonlinearities (mNEPv). The mNEPv is encountered in applications such as the computation of joint numerical radius of matrices, best rank-one approximation of third-order partial symmetric tensors, and distance to singularity for dissipative Hamiltonian differential-algebraic equations. We first present a variational characterization of the mNEPv. Based on the variational characterization, we provide a geometric interpretation of the self-consistent-field (SCF) iterations for solving the mNEPv, prove the global convergence of the SCF, and devise an accelerated SCF. Numerical examples from a variety of applications demonstrate the theoretical properties and computational efficiency of the SCF and its acceleration.

There are synergies of research interests and industrial efforts in modeling fairness and correcting algorithmic bias in machine learning. In this paper, we present a scalable algorithm for spectral clustering (SC) with group fairness constraints. Group fairness is also known as statistical parity where in each cluster, each protected group is represented with the same proportion as in the entirety. While FairSC algorithm (Kleindessner et al., 2019) is able to find the fairer clustering, it is compromised by high costs due to the kernels of computing nullspaces and the square roots of dense matrices explicitly. We present a new formulation of underlying spectral computation by incorporating nullspace projection and Hotelling's deflation such that the resulting algorithm, called s-FairSC, only involves the sparse matrix-vector products and is able to fully exploit the sparsity of the fair SC model. The experimental results on the modified stochastic block model demonstrate that s-FairSC is comparable with FairSC in recovering fair clustering. Meanwhile, it is sped up by a factor of 12 for moderate model sizes. s-FairSC is further demonstrated to be scalable in the sense that the computational costs of s-FairSC only increase marginally compared to the SC without fairness constraints.

We present a comprehensive convergence analysis for the self-consistent field (SCF) iteration to solve a class of nonlinear eigenvalue problems with eigenvector dependency (NEPvs). Using the tangent-angle matrix as an intermediate measure for approximation error, we establish new formulas for two fundamental quantities that characterize the local convergence behavior of the plain SCF: the local contraction factor and the local asymptotic average contraction factor. In comparison with previously established results, new convergence rate estimates provide much sharper bounds on the convergence speed. As an application, we extend the convergence analysis to a popular SCF variant---the level-shifted SCF. The effectiveness of the convergence rate estimates is demonstrated numerically for NEPvs arising from solving the Kohn--Sham equation in electronic structure calculation and the Gross--Pitaevskii equation for modeling of the Bose--Einstein condensation.

A thorough backward stability analysis of Hotelling's deflation, an explicit external deflation procedure through low-rank updates for computing many eigenpairs of a symmetric matrix, is presented. Computable upper bounds of the loss of the orthogonality of the computed eigenvectors and the symmetric backward error norm of the computed eigenpairs are derived. Sufficient conditions for the backward stability of the explicit external deflation procedure are revealed. Based on these theoretical results, the strategy for achieving numerical backward stability by dynamically selecting the shifts is proposed. Numerical results are presented to corroborate the theoretical analysis and to demonstrate the stability of the procedure for computing many eigenpairs of large symmetric matrices arising from applications.

We consider the optimization problem in which a continuous convex function is to be minimized over the joint numerical range of two Hermitian matrices. When those matrices are of large size, solving such problems by convex optimization can be computationally expensive. The goal of this paper is to present a novel nonlinear eigenvector method to accelerate the computation. We will show that the global minimizer of the optimization problem corresponds to a solution of a nonlinear eigenvalue problem with eigenvector nonlinearity (NEPv). The special structure of this NEPv allows for an efficient sequential subspace search algorithm, which is a nonlinear analogue to the NEPv of the commonly applied locally optimal conjugate gradient descent methods for Hermitian linear eigenvalue problems. Our new algorithm can be proven globally convergent to an eigenvector of the NEPv. Implementation details such as block iteration and preconditioning will be discussed. Numerical examples, with applications in computing the coercivity constant of boundary integral operators and solving multicast beamforming problems, show the effectiveness of our approach.

Some scientific and engineering applications need to compute a large number of eigenpairs of a large Hermitian matrix. Though the Lanczos method is effective for computing a few eigenvalues, it can be expensive for computing a large number of eigenpairs (e.g., in terms of computation and communication). To improve the performance of the method, in this paper, we study an s-step variant of thick-restart Lanczos (TRLan) combined with an explicit external deflation (EED). The s-step method generates a set of s basis vectors at a time and reduces the communication costs of generating the basis vectors. We then design a specialized matrix powers kernel (MPK) that reduces both the communication and computational costs by taking advantage of the special properties of the deflation matrix. We conducted numerical experiments of the new TRLan eigensolver using synthetic matrices and matrices from electronic structure calculations. The performance results on the Cori supercomputer at the National Energy Research Scientific Computing Center (NERSC) demonstrate the potential of the specialized MPK to significantly reduce the execution time of the TRLan eigensolver. The speedups of up to 3.1× and 5.3× were obtained in our sequential and parallel runs, respectively.

This paper presents a novel algorithm to compute the real stability radius for a linear delay system of retarded type with multiple delays. The real stability radius is the distance to instability, measured as the minimal real-valued perturbation that renders the system unstable. Our method is based on characterizing this distance to instability as the inverse of the global maximum of a real structured singular value function. We develop a criss-cross type algorithm that globally converges to this maximum, and whose convergence rate seems to be superlinear and sometimes quadratic in numerical experiments. The algorithm exploits that the intersections of these singular value functions with constant functions can be written as purely imaginary eigenvalues of certain delay eigenvalue problems (DEP) with positive and negative delays. This is an extension of the well-known linear case (without delays) where this results in algebraic eigenvalue problems. In addition, a novel numerical solver to compute all the imaginary eigenvalues of this DEP is also presented. It combines an approximation using spectral discretizations and an automatic procedure to determine the required number of discretization points. Finally, due to the presence of multiple eigenvalues at the maximum, these approximations are corrected with a block-Newton algorithm for nonlinear eigenvalue problems.

This paper is concerned with subspace acceleration techniques for computing the Crawford number, that is, the distance between zero and the numerical range of a matrix $A$. Our approach is based on an eigenvalue optimization characterization of the Crawford number. We establish local convergence of order $1+\sqrt{2}\approx 2.4$ for an existing subspace method applied to such and other eigenvalue optimization problems involving a Hermitian matrix that depends analytically on one parameter. For the particular case of the Crawford number, we show that the relevant part of the objective function is strongly concave. In turn, this enables us to develop a subspace method that only uses three-dimensional subspaces but still achieves global convergence and a local convergence that is at least quadratic. A number of numerical experiments confirm our theoretical results and reveal that the established convergence orders appear to be tight.

We study the robust Rayleigh quotient optimization problem where the data matrices of the Rayleigh quotient are subject to uncertainties. We propose to solve such a problem by exploiting its characterization as a nonlinear eigenvalue problem with eigenvector nonlinearity (NEPv). For solving the NEPv, we show that a commonly used iterative method can be divergent due to a wrong ordering of the eigenvalues. Two strategies are introduced to address this issue: a spectral transformation based on nonlinear shifting and a reformulation using second-order derivatives. Numerical experiments for applications in robust generalized eigenvalue classification, robust common spatial pattern analysis, and robust linear discriminant analysis demonstrate the effectiveness of the proposed approaches.

The real $\varepsilon$-pseudospectrum of a real matrix $A$ consists of the eigenvalues of all real matrices that are $\varepsilon$-close in spectral norm to $A$. The closeness is commonly measured in spectral or Frobenius norm. The real pseudospectral $\varepsilon$-pseudospectral abscissa, which is the largest real part of these eigenvalues for a prescribed value $\varepsilon$, measures the structured robust stability of $A$ w.r.t. real perturbations. In this paper, we introduce a criss-cross type algorithm to compute the real pseudospectral $\varepsilon$-pseudospectral abscissa for the spectral norm. Our algorithm is based on a superset characterization of the real pseudospectrum where each criss and cross search involves solving linear eigenvalue problems and singular value optimization problems. The new algorithm is proved to be globally convergent, and observed to be locally linearly convergent. Moreover, we propose a subspace projection framework in which we combine the criss-cross algorithm with subspace projection techniques to solve large-scale problems. The subspace acceleration is proved to be locally superlinearly convergent. The robustness and efficiency of the proposed algorithms are demonstrated on numerical examples.

The second-order Arnoldi (SOAR) procedure is an algorithm for computing an orthonormal basis of the second-order Krylov subspace. It has found applications in solving quadratic eigenvalue problems and model order reduction of second-order dynamical systems among others. Unfortunately, the SOAR procedure can be numerically unstable. The two-level orthogonal Arnoldi (TOAR) procedure has been proposed as an alternative to SOAR to cure the numerical instability. In this paper, we provide a rigorous stability analysis of the TOAR procedure. We prove that under mild assumptions, the TOAR procedure is backward stable in computing an orthonormal basis of the associated linear Krylov subspace. The benefit of the backward stability of TOAR is demonstrated by its high accuracy in structure-preserving model order reduction of second-order dynamical systems.

SummaryThe low‐rank damping term appears commonly in quadratic eigenvalue problems arising from physical simulations. To exploit the low‐rank damping property, we propose a Padé approximate linearization (PAL) algorithm. The advantage of the PAL algorithm is that the dimension of the resulting linear eigenvalue problem is only n + ℓm, which is generally substantially smaller than the dimension 2n of the linear eigenvalue problem produced by a direct linearization approach, where n is the dimension of the quadratic eigenvalue problem, and ℓ and m are the rank of the damping matrix and the order of a Padé approximant, respectively. Numerical examples show that by exploiting the low‐rank damping property, the PAL algorithm runs 33–47% faster than the direct linearization approach for solving modest size quadratic eigenvalue problems. Copyright © 2015 John Wiley & Sons, Ltd.