Abstract. The joint bidiagonalization (JBD) process of a regular matrix pair [Formula: see text] is mathematically equivalent to two joint Lanczos bidiagonalization processes of the upper and lower parts of the Q-factor of QR factorization of the stacked matrix [Formula: see text] when their starting vectors are closely related in a specific way. The resulting JBD method for computing extreme generalized singular values and corresponding left and right generalized singular vectors of [Formula: see text] realizes the standard Rayleigh–Ritz projection of the generalized singular value decomposition (GSVD) problem of [Formula: see text] onto the two left and one right subspaces generated by the JBD process. In this paper, the implicit restarting technique is nontrivially extended to the JBD process, and an implicitly restarted JBD (IRJBD) algorithm is developed with the shifts proposed and a few key implementation details addressed in finite precision arithmetic. Compact upper bounds are established for the residual norm of an approximate GSVD component in both exact and finite precision arithmetic, which are used to design efficient and reliable stopping criteria without explicitly computing approximate left and right generalized singular vectors at each iteration before convergence. Numerical experiments illustrate that IRJBD performs well and is at least competitive with the thick-restart JBD algorithm in terms of restarts and robustness. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SIMAX and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/jiazhongxiao/JBDmethod . [Formula: see text]
By assigning a nonnegative cost matrix to transitions between states, we propose an extension of the usual mean first passage times for a finite irreducible Markov chain to mean first passage costs. We then prove the existence of a (nonzero) fixed-point cost: i.e., a cost matrix which is, up to an overall multiplicative constant, the mean first passage cost matrix of itself. Moreover, we prove the uniqueness of the fixed-point cost, up to an overall multiplicative constant (and with the minor exception of an n-cycle permutation matrix). The entries of the fixed-point cost matrix thus define natural diffusion distances between the states since these distances are essentially invariant under taking expectations over random walks. We further show that the fixed-point cost matrix can be efficiently computed. In the last part of the paper, we return to general mean first passage costs and decompose them intro three components: one symmetric and two antisymmetric. One of the antisymmetric component vanishes in the case of the heat kernel-based diffusion transition matrix used by R. Coifman and S. Lafon. In this case, the symmetric component has an explicit connection to the Coifman--Lafon diffusion distance and isometric embedding.
Spectral methods are among the most widely used tools for analyzing data derived from real-world systems. Currently, many of the largest available datasets contain on the order of a billion or more vertices---and this number is only expected to grow. As a result, a central challenge is developing spectral methods, particularly eigenvalue solvers, that can handle the massive size and computational demands required to analyze such datasets. A notable characteristic of real-world datasets is the frequent presence of symmetries in the data's underlying graph structures. Our main theoretical contribution is to show that any graph (matrix)---whether weighted or unweighted, directed, or undirected---can be decomposed into a collection of smaller graphs (matrices), each defined with respect to the original graph's symmetries, and that this decomposition can be done while preserving the graph's (matrix's) full set of eigenvalues. As a direct application of this theory, we introduce the LEParD algorithm, which efficiently computes the eigenvalues of graphs exhibiting nontrivial symmetry---common in real-world data. We evaluate the LEParD algorithm on 110 datasets of varying types and sizes and find that, in approximately 82\% of the datasets, it outperforms standard eigenvalue solvers.
Orthogonal nonnegative matrix factorization (ONMF) extends traditional matrix factorization by incorporating nonnegativity and orthogonality constraints. It has gained significant attention for applications in dimension reduction and data clustering. This paper proposes a novel algorithm for solving ONMF within the block coordinate descent (BCD) framework. Both subsequence and whole sequence convergence of the algorithm are rigorously established, representing a major advancement in ONMF research. Extensive experiments on real-world datasets demonstrate the method's effectiveness and efficiency. Compared with several state-of-the-art ONMF techniques, the proposed algorithm achieves superior performance in terms of computational efficiency, relative error, and optimality gap.
Abstract. Dynamic mode decomposition (DMD), a fundamental methodology for data-driven dynamical systems analysis, faces three persistent limitations: restricted applicability to nonlinear systems, sensitivity to noise-induced instabilities, and computational inefficiency at scale. Recent advances partially address these constraints. However, these approaches still require a trade-off between computational efficiency and accuracy. This work introduces a unified randomized higher-order extended DMD (randomized HOEDMD) framework integrating randomized linear algebra with structured total least squares. Its innovations include Cholesky decomposition-enhanced randomized QB algorithms that reduce spatial complexity. Furthermore, we establish theoretical error bounds demonstrating quantifiable convergence to deterministic HOEDMD solutions under mild conditions. Evaluations across synthetic, cylinder wake flow, and functional magnetic resonance imaging (fMRI) datasets demonstrate (1) substantially improved computational efficiency while maintaining accuracy relative to the deterministic counterpart, and (2) superior precision compared to standard DMD and its randomized variants. Quantitatively, on the cylinder wake flow benchmark our method delivers [Formula: see text]42[Formula: see text] speedup at matched spectral accuracy; on large-scale voxel-level fMRI it avoids out-of-memory and completes in [Formula: see text]2.6 s. Overall, the proposed method provides an efficient, theoretically grounded approach for analyzing noise-contaminated, multiscale dynamical systems, effectively bridging computational tractability and dynamical fidelity.
Abstract. Many computational problems involve solving a linear system of equations, although only a subset of the entries of the solution are needed. In inverse problems, where the goal is to estimate unknown parameters from indirect noisy observations, it is not uncommon that the forward model linking the observed variables to the unknowns depends on variables that are not of primary interest, often referred to as nuisance parameters. In this article, we consider linear problems and propose a novel projection technique to eliminate, or at least mitigate, the contribution of the nuisance parameters in the model. We refer to this approach as spotlight inversion, as it allows to focus on only the portion of primary interest of the unknown parameter vector, leaving the uninteresting part in the shadow. The viability of the approach is illustrated with two computed examples, one where it works as model reduction for a finite element approximation of an elliptic PDE, the other amounting to local fanbeam X-ray tomography, spotlighting the region of interest that is part of the full target.
Abstract. Multigrid methods are powerful solvers for large-scale linear systems. Over the past decades, the convergence theory of multigrid methods for symmetric positive definite problems has been thoroughly studied and is currently well understood. However, fewer theoretical results exist for nonsymmetric problems. In [B. S. Southworth and T. A. Manteuffel, SIAM J. Matrix Anal. Appl., 45 (2024), pp. 1245–1258] the concept of so-called compatible transfer operators is introduced for nonsymmetric multigrid methods. Here we use this concept to establish a unified theoretical framework for multigrid methods. This framework contains different existing approaches not only for nonsymmetric methods but also for the symmetric case. We establish sharp convergence estimates, a criterion for optimal transfer (restriction and prolongation) operators is provided, and different convergence measures suitable for particular cases are explored. Moreover, our work generalizes previously well-known results.
The CANDECOMP/PARAFAC (CP) decomposition is widely used for analyzing multidimensional data, and the alternating least squares (CP-ALS) algorithm is a common method for its computation. CP rounding is the problem of computing a lower-rank CP decomposition of an input already in a higher-rank CP format. While the normal equations (NE) approach in CP-ALS is efficient for the CP rounding problem and frequently used, it becomes unstable in the presence of ill-conditioned subproblems. This paper presents a new QR-based CP-ALS method for CP rounding that preserves both numerical stability and computational efficiency. Our experiments show that the proposed method offers significant speedup over a previous QR-based approach and the Tensor Toolbox's NE-based implementation, particularly for higher-order tensors. Furthermore, our approach demonstrates a marked reduction in error for ill-conditioned problems, with error reductions several orders of magnitude smaller compared to the NE-based method, while achieving faster convergence and more accurate solutions. By using a more numerically stable approach, we can solve more problems in reduced working precision, which enables further reduction in time to solution.
In this paper, we originally derive the explicit expression for the optimal backward error of the multidimensional total least squares (mTLS) problem, with both full-rank and rankdeficient computed solutions considered. We propose the necessary and sufficient conditions for mTLS solutions, constructively characterizing the feasible set of the targeted perturbations. These the inherent difference in their problem structures. Theoretical analysis is still computationally prohibitive for large-scale applications. To address this issue, we propose computable estimates with tight bounds and further introduce randomized algorithms based on sketching techniques. Numerical experiments illustrate the effectiveness and the efficiency of the algorithms, offering a practical tool for iterative refinement. This work facilitates the backward stability analysis for the TLS problems.
We study conjugate partial-symmetric (CPS) tensors in this paper. The basis and dimension of the CPS tensor space is discussed and the conjugate tensor product decomposition of a CPS tensor is obtained from the CPS identifying matrix. We also introduce the conjugate partial-symmetric positive semidefinite (CPSPSD) tensors and the conjugate partial-symmetric completely positive (CPSCP) tensors in the CPS tensor space. The properties of the cone of CPSPSD tensors are given and the CPSB-eigenvalues and sum of squares representations of a CPSPSD tensor are investigated. In addition, we propose a semidefinite relaxation method to compute the smallest CPSB-eigenvalue. Further, we show that a CPSCP tensor has a CPSCP decomposition, where all real and imaginary parts of the decomposition vectors are nonnegative. A moment sum of squares (moment-SOS) relaxation algorithm is proposed for detecting whether a CPS tensor is CPSCP or not. If it is, we give a CPSCP decomposition for it; if it is not, we obtain a certificate for that. Low-rank CPSCP decompositions are also discussed.
The classical Birkhoff--von Neumann (BvN) decomposition expresses a given doubly stochastic matrix as a convex combination of permutation matrices. We investigate the BvN decomposition of symmetric doubly stochastic matrices where the permutation matrices in the decomposition are also symmetric, called symBvN decomposition. This decomposition is not always possible. Two pioneering theoretical works [M. W. Padb erg and L. A. Wolsey, Math. Programming, 29 (1984), pp. 1--14; V. V. Vazirani, An Extension of the Birkhoff--von Neumann Theorem to Non-Bipartite Graphs, preprint, arXiv:2010.05984, 2020] establish the conditions under which such a decomposition is possible using graph terminology. We propose a practical algorithm by combining these two works. A simple transformation converts any given symmetric doubly stochastic matrix, with possibly nonzero diagonal elements, to be the adjacency matrix of an edge-weighted undirected graph. The adjacency matrix of the resulting graph admits a symBvN decomposition if and only if the given matrix does so. The practicality of the proposed algorithm allows us to implement it, release its source code, and report the first set of experiments ever performed for the symBvN decomposition. Our experiments suggest that the proposed algorithm is as effective as the state-of-the-art algorithms for the classical BvN decomposition.
This paper introduces a deflation method for accurately computing the restricted singular value decomposition (RSVD) of a matrix triplet using its Neville representations (NRs). Our method exactly extracts all nonregular restricted singular values (RSVs) and deflates the triplet into a regular form while maintaining numerical accuracy. This approach ensures that all nonzero RSVs are computed with high relative accuracy, and all zero RSVs are exactly identified. It also exactly determines the dimensions of the associated nullspaces and their intersections. We establish rigorous perturbation bounds, demonstrating that the RSVD is accurately governed by its NR elements. The derived relative error bounds guarantee the high accuracy for the computed RSVD. Numerical experiments validate the claimed high relative accuracy.
We study the performance of the linear consensus algorithm on strongly connected directed graphs using the linear quadratic (LQ) cost as a performance measure. In particular, we derive bounds on the LQ cost by leveraging effective resistance and reversiblization. Our results extend previous analyses-which were limited to reversible cases-to the nonreversible setting. To facilitate this generalization, we introduce novel concepts, termed the back-and-forth path and the pivot node, which serve as effective alternatives to traditional techniques that require reversibility. Moreover, we apply our approach to Cayley graphs and random geometric graphs to estimate the LQ cost without the reversibility assumption. The proposed approach provides a framework that can be adapted to other contexts where reversibility is typically assumed.
The kth graph power G(k) of a graph G = (V, E) is the graph whose vertex set is V and in which two distinct vertices are adjacent if and only if their distance in G is at most k. The k-independence number alpha(k)(G) and distance-k chromatic number Xk(G) are then defined as the independence number and the chromatic number of G(k), respectively. We present a theoretical framework in which a wide range of bounds for the distance-k independence and chromatic numbers can be easily obtained and optimized in terms of the eigenvalues of G and a degree-k polynomial. We demonstrate the power of this method to derive sharp eigenvalue bounds for the two graph parameters. Moreover, we also show that several existing algebraic bounds fall in the proposed framework. Our approach is based on a combination of semidefinite programming and polynomial methods with spectral techniques.
In this paper, we investigate the robustness of frames in finite-dimensional real Hilbert spaces and their connection with totally nonsingular (TNS) matrices, providing novel characterizations and deterministic constructions for frames resilient to erasures. We establish that the robustness of a frame to erasures is intimately linked to the sparsity properties of vectors in the kernel of associated matrices. We provide explicit constructions of deterministic TNS matrices and full spark frames using structured matrices: Vandermonde matrices with negative bases, Cauchy matrices, and matrices with entries derived from distinct prime numbers or their square roots. Furthermore, we characterize when generalized Cauchy matrices can be orthogonal, enabling the construction of deterministic full spark tight frames over the real field---a significant advantage over complex-dependent Vandermonde construction.
Absolute value equations are nonlinear systems of the form Ax + B|x| = b, and they have seen rapid development in recent years. So far, the research has been focused primarily to the square case. In this paper, we aim to shift this focus and draw attention to overdetermined and underdetermined systems. For these rectangular cases, we address fundamental solvability issues. In particular, for overdetermined systems, we fully characterize cases where at most one solution exists. Similarly, for underdetermined systems, we present various conditions for the existence of at least one solution. The solvability of these systems is closely related to the full rank properties of certain sets of matrices. That is why we also investigate the full rank problems and derive different characterization conditions. Since absolute value equations are challenging to handle and many related problems are NP-hard, we examine rank-one (or more generally, fixed-rank) cases and demonstrate that certain questions can be answered efficiently. Finally, we pose several open problems, which are both easy to formulate and fundamental to solvability analysis.
The nonlinear GMRES (NGMRES) proposed by Washio and Oosterlee [Electron. Trans. Numer. Anal., 6 (1997), pp. 271--290] is an acceleration method for fixed point iterations. It has been demonstrated to be effective, but its convergence properties have not been extensively studied in the literature so far. In this work we aim to close some of this gap, by offering a convergence analysis for NGMRES applied to linear systems. A central part of our analysis focuses on identifying equivalences between NGMRES and the classical Krylov subspace GMRES method.
In this work, we propose a novel diagonalization-based preconditioner for the all-at-once linear system arising from the optimal control problem of parabolic equations. The proposed preconditioner is constructed based on an e-circulant modification to the rotated block diagonal (RBD) preconditioning technique and can be efficiently diagonalized by fast Fourier transforms in a parallel-in-time (PinT) fashion. To our knowledge, this marks the first application of the e-circulant modification to RBD preconditioning. Before our work, the studies of parallel-in-time preconditioning techniques for the optimal control problem are mainly focused on e-circulant modification to Schur complement based preconditioners, which involves multiplication of forward and backward evolutionary processes and thus square the condition number. Compared with those Schur complement based preconditioning techniques in the literature, the advantage of the proposed e-circulant modified RBD preconditioning is that it does not involve the multiplication of forward and backward evolutionary processes. The generalized minimal residual (GMRES) method is deployed on the preconditioned system, and we prove that the convergence rate of the preconditioned GMRES solver is independent of the matrix size and the regularization parameter when choosing e = \scrO(\surdT) with T being the temporal step-size. Numerical results are provided to demonstrate the effectiveness of our proposed solvers.
Kernel matrices are a key quantity in kernel-based approximation, and important properties such as stability and algorithmic convergence can be analyzed with their help. In this work, we refine a multivariate Ingham-type theorem, which is then leveraged to obtain novel and refined stability estimates on kernel matrices. For this, we focus on the case of finitely smooth kernels, such as the family of Mate'\rn or Wendland kernels, while noting that the results also extend to norm-equivalent kernels. In particular, we obtain results that relate the Rayleigh quotients of kernel matrices for kernels of different smoothness to each other. Finally, we comment on conclusions for the eigenvectors of these kernel matrices.
The aim of this work is to propose modified Wasserstein barycenters for probability measures defined on cartesian product sets which satisfy given marginal constraints. We focus on the specific case of Gaussian and Gaussian mixture distributions, as the proposed approach strongly relies on new results about properties of geometric means of covariance matrices. In the case of Gaussian distributions, the marginal-constrained modified Wasserstein barycenters can be analytically computed, while for Gaussian mixtures, computing the marginal-constrained barycenter consists in a postprocessing of the Gaussian mixture Wasserstein barycenter. In both cases, we provide numerical simulations illustrating the difference between Wasserstein barycenters and modified marginal-constrained Wasserstein barycenters. We moreover provide several test cases where the marginal-constrained Wasserstein barycenters interpolate better than regular Wasserstein barycenters, showcasing the practical interest of the proposed approach. As a by-product, we prove new results concerning marginal-preserving Wasserstein barycenters. Indeed, Wasserstein barycenters do not preserve marginals in general. In this work, as a consequence of the derived properties on the geometric mean of covariance matrices, we obtain sufficient and necessary conditions for the Wasserstein barycenter between two Gaussian distributions to preserve marginals and provide necessary conditions in the case of more than two Gaussians.
