Search

The GMRES algorithm of Saad and Schultz [SIAM J. Sci. Stat. Comput., 7 (1986), pp. 856-869] is an iterative method for approximately solving linear systems Ax = b, with initial guess x0 and residual r0 = b Ax0. The algorithm employs the Arnoldi process to generate the Krylov basis vectors (the columns of Vκ). It is well known that this process can be viewed as a QR factorization of the matrix Bκ = [ r0, AVκ] at each iteration. Despite an O (ε)κ (Bκ) loss of orthogonality, for unit roundoff ε and condition number κ, the modified Gram-Schmidt formulation was shown to be backward stable in the seminal paper by Paige et al. [SIAM J. Matrix Anal. Appl., 28 (2006), pp. 264-284]. We present an iterated Gauss-Seidel formulation of the GMRES algorithm (IGS-GMRES) based on the ideas of Ruhe [Linear Algebra Appl., 52 (1983), pp. 591-601] and Świrydowicz et al. [Numer. Linear Algebra Appl., 28 (2020), pp. 1-20]. IGS-GMRES maintains orthogonality to the level O (ε)κ (Bκ) or O (ε), depending on the choice of one or two iterations; for two Gauss-Seidel iterations, the computed Krylov basis vectors remain orthogonal to working accuracy and the smallest singular value of Vκ remains close to one. The resulting GMRES method is thus backward stable. We show that IGS-GMRES can be implemented with only a single synchronization point per iteration, making it relevant to large-scale parallel computing environments. We also demonstrate that, unlike MGS-GMRES, in IGS-GMRES the relative Arnoldi residual corresponding to the computed approximate solution no longer stagnates above machine precision even for highly nonnormal systems. [ABSTRACT FROM AUTHOR]

We consider the iterative method LSQR for solving minx ‖Ax-b‖2. LSQR is based on the Golub-Kahan bidiagonalization process and at every step produces an iterate that minimizes the norm of the residual vector over a Krylov subspace Kk. The 2-norm of the error is known to decrease monotonically, although it is not minimized over Kk. Given a lower bound on the smallest singular value of A, we show that in exact arithmetic the solution lies in the interior of a certain ellipsoid and that the LSQR iterate lies on the boundary of this ellipsoid. We use this result to derive new 2-norm error bounds for LSQR. Although our bounds are not much smaller than the existing ones, we show that they are sharp in the following sense: if the only information we use is our lower bound on σmin(A) plus the information gained by running k steps of LSQR, then our bounds cannot be improved. We also show how to choose a point with an error bound smaller than our corresponding bound for the LSQR error, although its true error is not necessarily smaller than the true LSQR error. As LSQR is formally equivalent to the conjugate gradient (CG) method applied to the normal equations ATAx = ATb, we derive analogous error bounds for CG. Our bounds for CG apply to any system Ax = b where A is symmetric positive definite. [ABSTRACT FROM AUTHOR]

We prove sharp lower bounds for the smallest singular value of a partial Fourier matrix with arbitrary "off the grid"" nodes (equivalently, a rectangular Vandermonde matrix with the nodes on the unit circle) in the case when some of the nodes are separated by less than the inverse bandwidth. The bound is polynomial in the reciprocal of the so-called superresolution factor, while the exponent is controlled by the maximal number of nodes which are clustered together. As a corollary, we obtain sharp minimax bounds for the problem of sparse superresolution on a grid under the partial clustering assumptions. [ABSTRACT FROM AUTHOR]

A sketch-and-select Arnoldi process to generate a well-conditioned basis of a Krylov space at low cost is proposed. At each iteration the procedure utilizes randomized sketching to select a limited number of previously computed basis vectors to project out of the current basis vector. The computational cost grows linearly with the dimension of the Krylov space. The subset selection problem for the projection step is approximately solved with a number of heuristic algorithms and greedy methods used in statistical learning and compressive sensing. [ABSTRACT FROM AUTHOR]

While projection-based reduced order models can reduce the dimension of full order solutions, the resulting reduced models may still contain terms that scale with the full order dimension. Hyper-reduction techniques are sampling-based methods that further reduce this computational complexity by approximating such terms with a much smaller dimension. The goal of this work is to introduce the points selection algorithm developed by Shin and Xiu [SIAM J. Sci. Comput., 38 (2016), pp. A385--A411] as a hyper-reduction method. The selection algorithm was originally proposed as a stochastic collocation method for uncertainty quantification. Since the algorithm aims at maximizing a quantity S that measures both the column orthogonality and the determinant, we refer to the algorithm as S-OPT. Numerical examples are provided to demonstrate the performance of S-OPT and to compare its performance with a gappy proper orthogonal decomposition (POD) algorithm. We found that using the S-OPT algorithm is shown to predict the full order solutions with higher accuracy than gappy POD especially when the number of sampling points is small, although we note that S-OPT shows slow asymptotic convergence with respect to the number of samples for some applications, e.g., Lagrangian hydrodynamics. [ABSTRACT FROM AUTHOR]

We propose a novel variant of the Localized Orthogonal Decomposition (LOD) method for addressing time-harmonic scattering problems of Helmholtz type with high wavenumber k. This method operates on a coarse mesh of width HH and identifies local finite element source terms that produce rapidly decaying responses under the solution operator. These source terms can be constructed with high accuracy from independent local snapshot solutions on patches of width ℓH, and they are used as problem-adapted basis functions. Compared to classical LOD and other state-of-the-art multiscale methods, our approach demonstrates that the localization error decays super-exponentially as the oversampling parameter ℓ increases. This indicates that optimal convergence is achieved under a substantially relaxed oversampling condition of ℓ≳(log⁡κ/H)(d−1)/d, where d represents the spatial dimension. Numerical experiments highlight the significant improvements in both offline and online performance of the method, even in the presence of heterogeneous media and perfectly matched layers. [ABSTRACT FROM AUTHOR]

This paper concerns the asymmetric transport observed along interfaces separating two-dimensional bulk topological insulators modeled by (continuous) differential Hamiltonians and how such asymmetry persists after numerical discretization. We first demonstrate that a relevant edge current observable is quantized and robust to perturbations for a large class of elliptic Hamiltonians. We then establish a bulk edge correspondence stating that the observable equals an integer-valued bulk difference invariant depending solely on the bulk phases. We next show how to extend such results to periodized Hamiltonians amenable to standard numerical discretizations. A form of no-go theorem implies that the asymmetric transport of periodized Hamiltonians necessarily vanishes. We introduce a filtered version of the edge current observable and show that it is approximately stable against perturbations and converges to its quantized limit as the size of the computational domain increases. To illustrate the theoretical results, we finally present numerical simulations that approximate the infinite domain edge current with high accuracy and show that it is approximately quantized even in the presence of perturbations. [ABSTRACT FROM AUTHOR]

We describe LNLQ for solving the least-norm problem min ||x|| subject to Ax = b, using the Golub-Kahan bidiagonalization of [b A]. Craig's method is known to be equivalent to applying the conjugate gradient method to the normal equations of the second kind (AAT y = b, x = AT y); LNLQ is equivalent to applying SYMMLQ. If an underestimate of the smallest singular value is available, error upper bounds for both x and y are available cheaply at each iteration. LNLQ is a companion method to the least-squares solver LSLQ [R. Estrin, D. Orban, and M. A. Saunders, SIAM J. Matrix Anal. Appl., 40 (2019b), pp. 235{253] which is equivalent to SYMMLQ on the conventional normal equations. We show that the error bounds are tight and comparable to the bounds suggested by Arioli [SIAM J. Matrix Anal. Appl., 34 (2013), pp. 571-592] for CRAIG. A sliding window technique allows us to tighten the error bound for y at the expense of a few additional scalar operations per iteration. We illustrate the tightness of the error bounds on two standard test problems and on the computation of an inexact gradient in the context of a penalty method for PDE-constrained optimization. [ABSTRACT FROM AUTHOR]

Reliable estimates for the condition number of a large, sparse, real matrix A are important in many applications. To get an approximation for the condition number κ(A), an approximation for the smallest singular value is needed. Standard Krylov subspaces are usually unsuitable for finding a good approximation to the smallest singular value. Therefore, we study extended Krylov subspaces which turn out to be ideal for the simultaneous approximation of both the smallest and largest singular value of a matrix. First, we develop a new extended Lanczos bidiagonalization method. With this method we obtain a lower bound for the condition number. Moreover, the method also yields probabilistic upper bounds for κ(A). The user can select the probability with which the upper bound holds, as well as the ratio of the probabilistic upper bound and the lower bound. [ABSTRACT FROM AUTHOR]

Whenever we use devices to take measurements, calibration is indispensable. While the purpose of calibration is to reduce bias and uncertainty in the measurements, it can be quite dificult, expensive, and sometimes even impossible to implement. We study a challenging problem called self-calibration, i.e., the task of designing an algorithm for devices so that the algorithm is able to perform calibration automatically. More precisely, we consider the setup y = A(d)x + ε where only partial information about the sensing matrix A(d) is known and where A(d) linearly depends on d. The goal is to estimate the calibration parameter d (resolve the uncertainty in the sensing process) and the signal/ object of interest x simultaneously. For three different models of practical relevance, we show how such a bilinear inverse problem, including blind deconvolution as an important example, can be solved via a simple linear least squares approach. As a consequence, the proposed algorithms are numerically extremely eficient, thus potentially allowing for real-time deployment. We also present a variation of the least squares approach, which leads to a spectral method, where the solution to the bilinear inverse problem can be found by computing the singular vector associated with the smallest singular value of a certain matrix derived from the bilinear system. Explicit theoretical guarantees and stability theory are derived for both techniques, and the number of sampling complexity is nearly optimal (up to a poly-log factor). Applications in imaging sciences and signal processing are discussed, and numerical simulations are presented to demonstrate the effectiveness and eficiency of our approach. [ABSTRACT FROM AUTHOR]

This paper describes a method that can be used either to compute exactly or toapproximate, in an appropriate norm, the distance of a matrix pair (A;B) from a pair ( ~ A; ~B) for which the product ~ A~B T has rank strictly less than the maximum of the ranks of A or B. If A and B both have less than full column rank, then the exact distance is given by the smallest singular value of an augmented matrix that is parameterized by a single real parameter. The particular value of the parameter that gives the distance can be determined by a line search. Even when A or B has full column rank, it is still always possible to at least approximate the distance. The method and its properties are demonstrated with several small numerical examples. [ABSTRACT FROM AUTHOR]

Any model order reduced dynamical system that evolves a modal decomposition to approximate the discretized solution of a stochastic PDE can be related to a vector field tangent to the manifold of fixed rank matrices. The dynamically orthogonal (DO) approximation is the canonical reduced-order model for which the corresponding vector field is the orthogonal projection of the original system dynamics onto the tangent spaces of this manifold. The embedded geometry of the fixed rank matrix manifold is thoroughly analyzed. The curvature of the manifold is characterized and related to the smallest singular value through the study of the Weingarten map. Differentiability results for the orthogonal projection onto embedded manifolds are reviewed and used to derive an explicit dynamical system for tracking the truncated singular value decomposition (SVD) of a time-dependent matrix. It is demonstrated that the error made by the DO approximation remains controlled under the minimal condition that the original solution stays close to the low rank manifold, which translates into an explicit dependence of this error on the gap between singular values. The DO approximation is also justified as the dynamical system that applies instantaneously the SVD truncation to optimally constrain the rank of the reduced solution. Riemannian matrix optimization is investigated in this extrinsic framework to provide algorithms that adaptively update the best low rank approximation of a smoothly varying matrix. The related gradient ow provides a dynamical system that converges to the truncated SVD of an input matrix for almost every initial datum. [ABSTRACT FROM AUTHOR]

Developing large-scale distributed methods that are robust to the presence of adversarial or corrupted workers is an important part of making such methods practical for real-world problems. In this paper, we propose an iterative approach that is adversary-tolerant for convex optimization problems. By leveraging simple statistics, our method ensures convergence and is capable of adapting to adversarial distributions. Additionally, the efficiency of the proposed methods for solving convex problems is shown in simulations with the presence of adversaries. Through simulations, we demonstrate the efficiency of our approach in the presence of adversaries and its ability to identify adversarial workers with high accuracy and tolerate varying levels of adversary rates. [ABSTRACT FROM AUTHOR]

We propose a new nonlinear preconditioned conjugate gradient (PCG) method in real arithmetic for computing the ground states of rotational Bose-Einstein condensate, modeled by the Gross--Pitaevskii equation. Our algorithm presents a few improvements of the PCG method in complex arithmetic studied by Antoine, Levitt, and Tang [J. Comput. Phys., 343 (2017), pp. 92-109]. We show that the special structure of the energy functional E(φ) and its gradient with respect to φ can be fully exploited in real arithmetic to evaluate them more efficiently. We propose a simple approach for fast evaluation of the energy functional, which enables exact line search. Most importantly, we derive the discrete Hessian operator of the energy functional and propose a shifted Hessian preconditioner for PCG, with which the ideal preconditioned Hessian has favorable eigenvalue distributions independent of the mesh size. This suggests that PCG with our ideal Hessian preconditioner is expected to exhibit mesh size-independent asymptomatic convergence behavior. In practice, our preconditioner is constructed by incomplete Cholesky factorization of the shifted discrete Hessian operator based on high-order finite difference discretizations. Numerical experiments in two-dimensional (2D) and three-dimensional (3D) domains show the efficiency of fast energy evaluation, the robustness of exact line search, and the improved convergence of PCG with our new preconditioner in iteration counts and runtime, notably for more challenging rotational BEC problems with high nonlinearity and rotational speed. [ABSTRACT FROM AUTHOR]

This work introduces a parallel and rank-adaptive matrix integrator for dynamical low-rank approximation. The method is related to the previously proposed rank-adaptive basis update and Galerkin (BUG) integrator but differs significantly in that all arising differential equations, both for the basis and the Galerkin coefficients, are solved in parallel. Moreover, this approach eliminates the need for a potentially costly coefficient update with augmented basis matrices. The integrator also incorporates a new step rejection strategy that enhances the robustness of both the parallel integrator and the BUG integrator. By construction, the parallel integrator inherits the robust error bound of the BUG and projector-splitting integrators. Comparisons of the parallel and BUG integrators are presented by a series of numerical experiments which demonstrate the efficiency of the proposed method, for problems from radiative transfer and radiation therapy. [ABSTRACT FROM AUTHOR]

In this work we propose, analyze, and test a new multiscale finite element method called Multiscale Hybrid (MH) method. The method is built as a close relative to the Multiscale Hybrid Mixed (MHM) method, but with the fundamental difference that a novel definition of the Lagrange multiplier is introduced. The practical implication of this is that both the local problems to compute the basis functions, as well as the global problem, are elliptic, as opposed to the MHM method (and also other previous methods) where a mixed global problem is solved and constrained local problems are solved to compute the local basis functions. The error analysis of the method is based on a hybrid formulation, and a static condensation process is done at the discrete level, so the final global system only involves the Lagrange multipliers. We tested the performance of the method by means of numerical experiments for problems with multiscale coefficients, and we carried out comparisons with the MHM method in terms of performance, accuracy, and memory requirements. [ABSTRACT FROM AUTHOR]

This paper is devoted to the development and convergence analysis of greedy reconstruction algorithms based on the strategy presented in [Y. Maday and J. Salomon, Joint Proceedings of the 48th IEEE Conference on Decision and Control and the 28th Chinese Control Conference, 2009, pp. 375--379]. These procedures allow the design of a sequence of control functions that ease the identification of unknown operators in nonlinear dynamical systems. The original strategy of greedy reconstruction algorithms is based on an offline/online decomposition of the reconstruction process and an ansatz for the unknown operator obtained by an a priori chosen set of linearly independent matrices. In the previous work [S. Buchwald, G. Ciaramella, and J. Salomon, SIAM J. Control Optim., 59 (2021), pp. 4511--4537], convergence results were obtained in the case of linear identification problems. We tackle here the more general case of nonlinear systems. More precisely, we introduce a new greedy algorithm based on the linearized system. We show that the controls obtained with this new algorithm lead to the local convergence of the classical Gauss--Newton method applied to the online nonlinear identification problem. We then extend this result to the controls obtained on nonlinear systems where a local convergence result is also proved. The main convergence results are obtained for dynamical systems with linear and bilinear control structures. [ABSTRACT FROM AUTHOR]

Classical primal-dual algorithms attempt to solve max∣1 miι⅛ £(æ, μ) by alternately minimizing over the primal variable x through primal descent and maximizing the dual variable μ through dual ascent. However, when Z2(x, μ) is highly nonconvex with complex constraints in x, the minimization over x may not achieve global optimality and, hence, the dual ascent step loses its valid intuition. This observation motivates us to propose a new class of primal-dual algorithms for nonconvex constrained optimization with the key feature to reverse dual ascent to a conceptually new dual descent, in a sense, elevating the dual variable to the same status as the primal variable. Surprisingly, this new dual scheme achieves some best iteration complexities for solving nonconvex optimization problems. In particular, when the dual descent step is scaled by a fractional constant, we name it scaled dual descent (SDD), otherwise, unsealed dual descent (UDD). For nonconvex multiblock optimization with nonlinear equality constraints, we propose SDD-alternating direction method of multipliers (SDD-ADMM) and show that it finds an e-stationary solution in O{e~4) iterations. The complexity is further improved to O{e~3) and O{e~2) under proper conditions. We also propose UDD-augmented Lagrangian method (UDD-ALM), combining UDD with ALM, for weakly convex minimization over affine constraints. We show that UDD-ALM finds an e-stationary solution in O{e~2) iterations. These complexity bounds for both algorithms either achieve or improve the best-known results in the ADMM and ALM literature. Moreover, SDD-ADMM addresses a longstanding limitation of existing ADMM frameworks. [ABSTRACT FROM AUTHOR]

Sketch-and-precondition techniques are efficient and popular for solving large least squares (LS) problems of the form Ax = b with A∈ Rm×n and m≫n. This is where A is "sketched" to a smaller matrix SA with S∈ R[cn]×m for some constant c > 1 before an itera- tive LS solver computes the solution to Ax = b with a right preconditioner P, where P is constructed from SA. Prominent sketch-and-precondition LS solvers are Blendenpik and LSRN. We show that the sketch-and-precondition technique in its most commonly used form is not numerically stable for ill-conditioned LS problems. For provable and practical backward stability and optimal residuals, we suggest using an unpreconditioned iterative LS solver on (AP)z = b with x = Pz. Provided the condition number of A is smaller than the reciprocal of the unit roundoff, we show that this modification ensures that the computed solution has a backward error comparable to the iterative LS solver applied to a well-conditioned matrix. Using smoothed analysis, we model floating-point rounding errors to argue that our modification is expected to compute a backward stable solution even for arbitrarily ill-conditioned LS problems. Additionally, we provide experimental evidence that using the sketch-and-solve solution as a starting vector in sketch-and-precondition algorithms (as suggested by Rokhlin and Tygert in 2008) should be highly preferred over the zero vector. The initialization often results in much more accurate solutions albeit not always backward stable ones. [ABSTRACT FROM AUTHOR]

We derive closed formulas for the condition number of a linear function of the total least-squares solution. Given an overdetermined linear system Ax ≈ b, we show that this condition number can be computed using the singular values and the right singular vectors of [A, b] and A. We also provide an upper botlnd that requires the computation of the largest and the smallest singular value of [A, b] and the smallest singular value of A. In numerical examples, we compare these values and the resulting forward error bounds with the error estimates given by Van Huffel and Vandewalle [The Total Least Squares Problem: Computational Aspects and Analysis, Frontiers Appl. Math. 9, SIAM, Philadelphia, 1991], and we show the limitation of the first order approach. [ABSTRACT FROM AUTHOR]

We provide expressions for the residual norms when using the full orthogonalization method and the generalized minimum residual method for solving linear systems. They involve a triangular submatrix of the Hessenberg matrix generated by the Arnoldi process. This allows one to obtain bounds showing that the norm of the residual decreases to zero when the smallest singular value of this triangular matrix goes to zero. Numerical examples show that even though these bounds are not sharp they describe quite well the rate of decrease of the residual norm. [ABSTRACT FROM AUTHOR]

We define in the space of n×m matrices of rank n, n ≤ m, the condition Riemannian structure as follows: For a given matrix A the tangent space at A is equipped with the Hermitian inner product obtained by multiplying the usual Frobenius inner product by the inverse of the square of the smallest singular value of A denoted σn(A). When this smallest singular value has multiplicity 1, the function A → log(σn(A)−2) is a convex function with respect to the condition Riemannian structure that is t → log(σn(A(t))−2) is convex, in the usual sense for any geodesic A(t). In a more abstract setting, a function α defined on a Riemannian manifold (M, ± ) is said to be self-convex when log α(γ(t)) is convex for any geodesic in (M, α ±, ). Necessary and sufficient conditions for self-convexity are given when α is C2. When α(x) = d(x,N)−2, where d(x,N) is the distance from x to a C2 submanifold N ⊂Rj, we prove that α is self-convex when restricted to the largest open set of points x where there is a unique closest point in N to x. We also show, using this more general notion, that the square of the condition number ∙A∙F /σn(A) is self-convex in projective space and the solution variety. [ABSTRACT FROM AUTHOR]

Randomized Kaczmarz is a simple iterative method for finding solutions of linear systems Ax=b. We point out that the arising sequence (xk)∞k=1 tends to converge to the solution xin an interesting way: generically, ask→∞, xk−x tends to the singular vector of Acorresponding to the smallest singular value. This has interesting consequences: in particular, the error analysis of Strohmer and Vershyn in is optimal. It also quantifies the “preconvergence” phenomenon where the method initially seems to converge faster. This fact also allows for a fast computation of vectors x for which the Rayleigh quotient‖Ax‖/‖x‖is small: solve Ax= 0 via randomized Kaczmarz. [ABSTRACT FROM AUTHOR]

Computational methods for thermal radiative transfer problems exhibit high computational costs and a prohibitive memory footprint when the spatial and directional domains are finely resolved. A strategy to reduce such computational costs is dynamical low-rank approximation (DLRA), which represents and evolves the solution on a low-rank manifold, thereby significantly decreasing computational and memory requirements. Efficient discretizations for the DLRA evolution equations need to be carefully constructed to guarantee stability while enabling mass conservation. In this work, we focus on the Su--Olson closure leading to a linearized internal energy model and derive a stable discretization through an implicit coupling of internal energy and particle density. Moreover, we propose a rank-adaptive strategy to preserve local mass conservation. Numerical results are presented which showcase the accuracy and efficiency of the proposed low-rank method compared to the solution of the full system. [ABSTRACT FROM AUTHOR]

The Nyström method is a convenient strategy to obtain low-rank approximations to kernel matrices in nearly linear complexity. Existing studies typically use the method to approximate positive semidefinite matrices with low or modest accuracies. In this work, we propose a series of heuristic strategies to make the Nyström method reach high accuracies for nonsymmetric and/or rectangular matrices. The resulting methods (called high-accuracy Nyström methods) treat the Nyström method and a skinny rank-revealing factorization as a fast pivoting strategy in a progressive alternating direction refinement process. Two refinement mechanisms are used: alternating the row and column pivoting starting from a small set of randomly chosen columns, and adaptively increasing the number of samples until a desired rank or accuracy is reached. A fast subset update strategy based on the progressive sampling of Schur complements is further proposed to accelerate the refinement process. Efficient randomized accuracy control is also provided. Relevant accuracy and singular value analysis is given to support some of the heuristics. Extensive tests with various kernel functions and data sets show how the methods can quickly reach prespecified high accuracies in practice, sometimes with quality close to SVDs, using only small numbers of progressive sampling steps. [ABSTRACT FROM AUTHOR]

Quadratic minimization problems with orthogonality constraints (QMPO) play an important role in many applications of science and engineering. However, some existing methods may suffer from low accuracy or heavy workload for large-scale QMPO. Krylov subspace methods are popular for large-scale optimization problems. In this work, we propose a block Lanczos method for solving the large-scale QMPO. In the proposed method, the original problem is projected into a small-sized one, and the Riemannian trust-region method is employed to solve the reduced QMPO. Convergence results on the optimal solution, the optimal objective function value, the multiplier, and the KKT error are established. Moreover, we give the convergence speed of the approximate solution and show that if the block Lanczos process terminates, then an exact KKT solution is derived. Numerical experiments illustrate the numerical behavior of the proposed algorithm and demonstrate that it is more powerful than many state-of-the-art algorithms for large-scale QMPO. [ABSTRACT FROM AUTHOR]

The Tucker tensor decomposition is a natural extension of the singular value decomposition (SVD) to multiway data. We propose to accelerate Tucker tensor decomposition algorithms by using randomization and parallelization. We present two algorithms that scale to large data and many processors, significantly reduce both computation and communication cost compared to previous deterministic and randomized approaches, and obtain nearly the same approximation errors. The key idea in our algorithms is to perform randomized sketches with Kronecker-structured random matrices, which reduces computation compared to unstructured matrices and can be implemented using a fundamental tensor computational kernel. We provide probabilistic error analysis of our algorithms and implement a new parallel algorithm for the structured randomized sketch. Our experimental results demonstrate that our combination of randomization and parallelization achieves accurate Tucker decompositions much faster than alternative approaches. We observe up to a 16X speedup over the fastest deterministic parallel implementation on 3D simulation data. [ABSTRACT FROM AUTHOR]

We develop two new sets of stable, rank-adaptive dynamically orthogonal Runge--Kutta (DORK) schemes that capture the high-order curvature of the nonlinear low-rank manifold. The DORK schemes asymptotically approximate the truncated singular value decomposition at a greatly reduced cost while preserving mode continuity using newly derived retractions. We show that arbitrarily high-order optimal perturbative retractions can be obtained, and we prove that these new retractions are stable. In addition, we demonstrate that repeatedly applying retractions yields a gradient-descent algorithm on the low-rank manifold that converges superlinearly when approximating a low-rank matrix. When approximating a higher-rank matrix, iterations converge linearly to the best low-rank approximation. We then develop a rank-adaptive retraction that is robust to overapproximation. Building off of these retractions, we derive two rank-adaptive integration schemes that dynamically update the subspace upon which the system dynamics are projected within each time step: the stable, optimal dynamically orthogonal Runge--Kutta (so-DORK) and gradient-descent dynamically orthogonal Runge--Kutta (gd-DORK) schemes. These integration schemes are numerically evaluated and compared on an ill-conditioned matrix differential equation, an advection-diffusion partial differential equation, and a nonlinear, stochastic reaction-diffusion partial differential equation. Results show a reduced error accumulation rate with the new stable, optimal and gradient-descent integrators. In addition, we find that rank adaptation allows for highly accurate solutions while preserving computational efficiency. [ABSTRACT FROM AUTHOR]

We present an online algorithm for time-varying semidefinite programs (TV-SDPs), based on the tracking of the solution trajectory of a low-rank matrix factorization, also known as the Burer--Monteiro factorization, in a path-following procedure. There, a predictor-corrector algorithm solves a sequence of linearized systems. This requires the introduction of a horizontal space constraint to ensure the local injectivity of the low-rank factorization. The method produces a sequence of approximate solutions for the original TV-SDP problem, for which we show that they stay close to the optimal solution path if properly initialized. Numerical experiments for a time-varying max-cut SDP relaxation demonstrate the computational advantages of the proposed method for tracking TV-SDPs in terms of runtime compared to off-the-shelf interior-point methods. [ABSTRACT FROM AUTHOR]

Anderson Acceleration (AA) is a popular acceleration technique to enhance the convergence of fixed-point iterations. The analysis of AA approaches typically focuses on the convergence behavior of a corresponding fixed-point residual, while the behavior of the underlying objective function values along the accelerated iterates is currently not well understood. In this paper, we investigate local properties of AA with restarting applied to a basic gradient scheme in terms of function values. Specifically, we show that AA with restarting is a local descent method and that it can decrease the objective function faster than the gradient method. These new results theoretically support the good numerical performance of AA when heuristic descent conditions are used for globalization and they provide a novel perspective on the convergence analysis of AA that is more amenable to nonconvex optimization problems. Numerical experiments are conducted to illustrate our theoretical findings. [ABSTRACT FROM AUTHOR]

We study a new two-time-scale stochastic gradient method for solving optimization problems, where the gradients are computed with the aid of an auxiliary variable under samples generated by time-varying Markov random processes controlled by the underlying optimization variable. These time-varying samples make gradient directions in our update biased and dependent, which can potentially lead to the divergence of the iterates. In our two-time-scale approach, one scale is to estimate the true gradient from these samples, which is then used to update the estimate of the optimal solution. While these two iterates are implemented simultaneously, the former is updated "faster" (using bigger step sizes) than the latter (using smaller step sizes). Our first contribution is to characterize the finite-time complexity of the proposed two-time-scale stochastic gradient method. In particular, we provide explicit formulas for the convergence rates of this method under different structural assumptions, namely, strong convexity, the Polyak--Lojasiewicz condition, and general nonconvexity. We apply our framework to policy optimization problems in control and reinforcement learning. First, we look at the infinite-horizon average-reward Markov decision process with finite state and action spaces and derive a convergence rate of Ő(k-2/5) for the online actor-critic algorithm under function approximation, which recovers the best known rate derived specifically for this problem. Second, we study the linear-quadratic regulator and show that an online actor-critic method converges with rate Ő(k-2/2). Third, we use the actor-critic algorithm to solve the policy optimization problem in an entropy regularized Markov decision process, where we also establish a convergence of Ő(k-2/3). The results we derive for both the second and third problems are novel and previously unknown in the literature. Finally, we briefly present the application of our framework to gradient-based policy evaluation algorithms in reinforcement learning. [ABSTRACT FROM AUTHOR]

Intensity-based image registration is critical for neuroimaging tasks, such as 3D reconstruction, times-series alignment, and common coordinate mapping. The gradient-based optimization methods commonly used to solve this problem require a careful selection of step-length. This limitation imposes substantial time and computational costs. Here we propose a gradient-independent rigid-motion registration algorithm based on the majorization-minimization (MM) principle. Each iteration of our intensity-based MM algorithm reduces to a simple point-set rigid registration problem with a closed form solution that avoids the step-length issue altogether. The details of the algorithm are presented, and an error bound for its more practical truncated form is derived. The performance of the MM algorithm is shown to be more effective than gradient descent on simulated images and Nissl stained coronal slices of mouse brain. We also compare and contrast the similarities and differences between the MM algorithm and another gradient-free registration algorithm called the block-matching method. Finally, extensions of this algorithm to more complex problems are discussed. [ABSTRACT FROM AUTHOR]

The cross-product matrix-based CJ-FEAST SVDsolver proposed previously by the authors is shown to compute the left singular vector possibly much less accurately than the right singular vector and may be numerically backward unstable when a desired singular value is small. In this paper, an alternative augmented matrix-based CJ-FEAST SVDsolver is proposed to compute the singular triplets of a large matrix A with the singular values in an interval [a, b] contained in the singular spectrum. The new CJ-FEAST SVDsolver is a subspace iteration applied to an approximate spectral projector of the augmented matrix [0, AT; A, 0] associated with the eigenvalues in [a, b], and it constructs approximate left and right singular subspaces independently, onto which A is projected to obtain the Ritz approximations to the desired singular triplets. Compact estimates are given for the accuracy of the approximate spectral projector constructed by the Chebyshev--Jackson series expansion in terms of series degree, and a number of convergence results are established. The new solver is proved to be always numerically backward stable. A convergence comparison of the cross-product-based and augmented matrix-based CJ-FEAST SVDsolvers is made, and a generalpurpose choice strategy between the two solvers is proposed for the robustness and overall efficiency. Numerical experiments confirm all the results and meanwhile demonstrate that the proposed solver is more robust and substantially more efficient than the corresponding contour integral-based versions that exploit the trapezoidal rule and the Gauss--Legendre quadrature to construct an approximate spectral projector. [ABSTRACT FROM AUTHOR]

Given a real matrix A with n columns, the problem is to approximate the Gram product AAT by c ≪ n weighted outer products of columns of A. Necessary and sufficient conditions for the exact computation of AAT (in exact arithmetic) from c = rank(A) columns depend on the right singular vector matrix of A. For a Monte Carlo matrix multiplication algorithm by Drineas et al. that samples outer products, we present probabilistic bounds for the two-norm relative error due to randomization. The bounds depend on the stable rank or the rank of A, but not on the matrix Given a real matrix A with n columns, the problem is to approximate the Gram product AAT by c ≥ n weighted outer products of columns of A. Necessary and sufficient conditions for the exact computation of AAT (in exact arithmetic) from c = rank(A) columns depend on the right singular vector matrix of A. For a Monte Carlo matrix multiplication algorithm by Drineas et al. that samples outer products, we present probabilistic bounds for the two-norm relative error due to randomization. The bounds depend on the stable rank or the rank of A, but not on the matrix dimensions. Numerical experiments illustrate that the bounds are informative, even for stringent success probabilities and matrices of small dimension. We also derive bounds for the smallest singular value and the condition number of matrices obtained by sampling rows from orthonormal matrices. [ABSTRACT FROM AUTHOR]

Canonical polyadic decomposition (CPD) is prevalent in chemometrics, signal processing, data mining, and many more fields. While many algorithms have been proposed to compute the CPD, alternating least squares (ALS) remains one of the most widely used algorithms for computing the decomposition. Recent works have introduced the notion of eigenvalues and singular values of a tensor and explored applications of eigenvectors and singular vectors in signal processing, data analytics, and various other fields. We introduce a new formulation for deriving singular values and vectors of a tensor by considering the critical points of a function differently from previous works. Computing these critical points in an alternating manner motivates an alternating optimization algorithm which corresponds to the ALS algorithm in the matrix case. However, for tensors with order greater than or equal to 3, it minimizes an objective function which is different from the commonly used least squares loss. Alternating optimization of this new objective leads to simple updates to the factor matrices with the same asymptotic computational cost as ALS. We show that a subsweep of this algorithm can achieve a superlinear convergence rate for exact CPD when the known rank is not larger than the mode lengths of the input tensor. We verify our theoretical arguments experimentally. We then view the algorithm as optimizing a Mahalanobis distance with respect to each factor with the ground metric dependent on the other factors. This perspective allows us to generalize our approach to interpolate between updates corresponding to the ALS and the new algorithm to manage the tradeoff between stability and fitness of the decomposition. Our experimental results show that for approximating synthetic and real-world tensors, this algorithm and its variants converge to a better conditioned decomposition with comparable and sometimes better fitness as compared to the ALS algorithm. [ABSTRACT FROM AUTHOR]

We consider the problem of reconstructing one-dimensional point sources from their Fourier measurements in a bounded interval [\Omega, \Omega]. This problem is known to be challenging in the regime where the spacing of the sources is below the Rayleigh length \pi \Omega. In this paper, we propose a superresolution algorithm, called iterative focusing-localization and filtering, to resolve closely spaced point sources from their multiple measurements that are obtained by using multiple unknown illumination patterns. The new proposed algorithm has a distinct feature in that it reconstructs the point sources one by one in an iterative manner and hence requires no prior information about the source numbers. The new feature also allows for a subsampling strategy that can reconstruct sources using small-sized Hankel matrices and thus circumvent the computation of singular-value decomposition for large matrices as in the usual subspace methods. In addition, the algorithm can be paralleled. A theoretical analysis of the methods behind the algorithm is also provided. The derived results imply a phase transition phenomenon in the reconstruction of source locations which is confirmed in the numerical experiment. Numerical results show that the algorithm can achieve a stable reconstruction for point sources with a minimum separation distance that is close to the theoretical limit. The efficiency and robustness of the algorithm have also been tested. This algorithm can be generalized to higher dimensions. [ABSTRACT FROM AUTHOR]

By introducing a new operator theory, we provide a unified mathematical theory for general source resolution in the multi-illumination imaging problem. Our main idea is to transform multiillumination imaging into single-snapshot imaging with a new imaging kernel that depends on both the illumination patterns and the point spread function of the imaging system. We therefore prove that the resolution of multi-illumination imaging is approximately determined by the essential cutoff frequency of the new imaging kernel, which is roughly limited by the sum of the cutoff frequency of the point spread function and the maximum essential frequency in the illumination patterns. Our theory provides a unified way to estimate the resolution of various existing super-resolution modalities and results in the same estimates as those obtained in experiments. In addition, based on the reformulation of the multi-illumination imaging problem, we also estimate the resolution limits for resolving both complex and positive sources by sparsity-based approaches. We show that the resolution of multi-illumination imaging is approximately determined by the new imaging kernel from our operator theory and better resolution can be realized by sparsity-promoting techniques in practice but only for resolving very sparse sources. This explains experimentally observed phenomena in some sparsity-based super-resolution modalities. [ABSTRACT FROM AUTHOR]

Convergence analysis of block iterative solvers for Hermitian eigenvalue problems and closely related research on properties of matrix-based signal filters are challenging and are attracting increased attention due to their recent applications in spectral data clustering and graph-based signal processing. We combine majorization-based techniques pioneered for investigating the Rayleigh--Ritz method in [A. V. Knyazev and M. E. Argentati, SIAM J. Matrix Anal. Appl., 31 (2010), pp. 1521-1537] with tools of classical analysis of the block power method by Rutishauser [Numer. Math., 13 (1969), pp. 4-13] to derive sharp convergence rate bounds of abstract block iterations, wherein tuples of tangents of principal angles or relative errors of Ritz values are bounded using majorization in terms of arranged partial sums and tuples of convergence factors. Our novel bounds are robust in the presence of clusters of eigenvalues, improve previous results, and are applicable to most known block iterative solvers and matrix-based filters, e.g., to block power, Chebyshev, and Lanczos methods combined with polynomial filtering. The sharpness of our bounds is fundamental, implying that the bounds cannot be improved without further assumptions. [ABSTRACT FROM AUTHOR]

A common approach for compressing large-scale data is through matrix sketching. In this work, we consider the problem of recovering low-rank matrices from two noisy linear sketches using the double sketching scheme discussed in Fazel et al. [Compressed sensing and robust recovery of low rank matrices, in Proceedings of the 42nd IEEE Asilomar Conference on Signals, Systems and Computers, 2008, pp. 1043-1047], which is based on an approach by Woolfe et al. [Appl. Comput. Harmon. Anal., 25 (2008), pp. 335-366]. Using tools from nonasymptotic random matrix theory, we provide the first theoretical guarantees characterizing the error between the output of the double sketch algorithm and the ground truth low-rank matrix. We apply our result to the problems of low-rank matrix approximation and low-tubal-rank tensor recovery. [ABSTRACT FROM AUTHOR]

We introduce a new method for constructing reduced-order models that accurately approximate the frequency domain behavior of large-scale, parametrized, linear dynamical systems. Using multivariate rational interpolation, we compute reduced-order models that match transfer function measurements of the large-scale system. The main tools are a new barycentric formula and Loewner matrices formed directly from measurements. This data-driven approach introduces a new degree of freedom for parametrized model reduction--the ability to choose separate reduced orders for each parameter. More precisely, each reduced order is determined by computing the rank of appropriate Loewner matrices. Moreover, we show how to control the pointwise approximation error through a new formula involving the barycentric form and the smallest singular value of a Loewner matrix. Finally, we also give new state-space realizations for multivariate rational functions and conclude with several numerical examples showcasing the effectiveness of the new method. [ABSTRACT FROM AUTHOR]

We propose a new interpretation of the generalized overdetermined eigenvalue problem (A -- λB) v ≈ 0 for two mxn (m > n) matrices A and B, its stability analysis, and an efficient algorithm for solving it. Usually, the matrix pencil {A-λB} does not have any rank deficient member. Therefore we aim to compute λ for which A - λB is as close as possible to rank deficient; i.e., we search for λ that locally minimize the smallest singular value over the matrix pencil {A - λB}. Practically, the proposed algorithm requires O(mn²) operations for computing all the eigenpairs. We also describe a method to compute practical starting eigenpairs. The effectiveness of the new approach is demonstrated with numerical experiments. A MATLAB-based implementation of the proposed algorithm can be found at http://www.mat.univie.ac.at/~neum/software/oeig/. [ABSTRACT FROM AUTHOR]

The solution of the total least squares (TLS) problems, minE,f ¦ ¦(E,f)¦ ¦F subject to (A + E)χ = b + f, can in the generic case be obtained from the right singular vector corresponding to the smallest singular value σn+1 of (A, b). When A is large and sparse (or structured) a method based on Rayleigh quotient iteration (RQI) has been suggested by Björck. In this method the problem is reduced to the solution of a sequence of symmetric, positive definite linear systems of the form (AT A - &isgma;¯²I)z = g, where &isgma¯ is an approximation to σn+1 . These linear systems are then solved by a preconditioned conjugate gradient method (PCGTLS). For TLS problems where A is large and sparse a (possibly incomplete) Cholesky factor of AT A can usually be computed, and this provides a very efficient preconditioner. The resulting method can be used to solve a much wider range of problems than it is possible to solve by using Lanczos-type algorithms directly for the singular value problem. In this paper the RQI-PCGTLS method is further developed, and the choice of initial approximation and termination criteria are discussed. Numerical results confirm that the given algorithm achieves rapid convergence and good accuracy. [ABSTRACT FROM AUTHOR]

Noise poses a challenge for learning dynamical-system models because already small variations can distort the dynamics described by trajectory data. This work builds on operator inference from scientific machine learning to infer low-dimensional models from high-dimensional state trajectories polluted with noise. The presented analysis shows that, under certain conditions, the inferred operators are unbiased estimators of the well-studied projection-based reduced operators from traditional model reduction. Furthermore, the connection between operator inference and projection-based model reduction enables bounding the mean-squared errors of predictions made with the learned models with respect to traditional reduced models. The analysis also motivates an active operator inference approach that judiciously samples high-dimensional trajectories with the aim of achieving a low mean-squared error by reducing the effect of noise. Numerical experiments with high-dimensional linear and nonlinear state dynamics demonstrate that predictions obtained with active operator inference have orders of magnitude lower mean-squared errors than operator inference with traditional, equidistantly sampled trajectory data. [ABSTRACT FROM AUTHOR]

In this paper we discuss a projection model order reduction method for a class of parametric linear evolution PDEs, which is based on the application of the Laplace transform. The main advantage of this approach is that the Laplace transform allows us to compute the solution directly at a given instant. This can be done efficiently by approximating the contour integral associated to the inverse Laplace transform by a suitable quadrature formula. This feature turns out to be very important in terms of classical model order reduction methodologies. In fact, differently from time stepping integrators (like Runge-Kutta methods), the use of the Laplace transform determines a significant reduction of the size of the reduced space. In this article we propose two reduction methods; the first one constructs a unique reduced space, based on the vectors representing the Laplace transform at the chosen quadrature points for different values of the parameters, and the second one instead constructs several reduced spaces, one at each quadrature node. This last approach is particularly well suited to a parallel computation approach, allowing a faster computation. We show the effectiveness of the method by some illustrative parabolic PDEs arising from finance. [ABSTRACT FROM AUTHOR]

A new approach to solving eigenvalue optimization problems for large structured matrices is proposed and studied. The class of optimization problems considered is related to computing structured pseudospectra and their extremal points, and to structured matrix nearness problems such as computing the distance to instability or to singularity under structured perturbations. The perturbation structure can be a general linear structure. We focus on the practically important cases of large matrices with a given sparsity pattern and on perturbation matrices with given range and corange. It is known that analogous eigenvalue optimization for unstructured complex matrices favorably works with rank-1 matrices. The novelty of the present paper is that structured eigenvalue optimization can still be performed with rank-1 matrices, which yields a significant reduction of storage and in some cases of the computational cost. Optimizers are shown to be rank-1 matrices orthogonally projected onto the given structure. This is used in the numerical algorithms. The numerical approach comes in two different versions. If the rank-1 matrices projected onto the structure form a manifold and if the orthogonal projection onto the corresponding tangent spaces is known and computationally inexpensive, as is the case for matrices with given range and corange, then the method relies on a gradient system on this manifold. Otherwise, in particular for matrices with a given sparsity pattern, the method relies on the orthogonal projection of the free gradient onto the tangent space of the manifold of complex rank-1 matrices. The solution of this rank-1 system is then projected onto the structure. It is shown that there is a bijective correspondence between the stationary points of the gradient system on the structure space and of the rank-1 system. Near a local minimizer the rank-1 tangent projection is very close to the identity map, and so the computationally favorable rank-1 projected system behaves locally like the gradient system. Numerical experiments illustrate the two approaches for large matrices with a given sparsity pattern and for perturbation matrices with given range and corange. [ABSTRACT FROM AUTHOR]

The problem of interest is the minimization of a nonlinear function subject to nonlinear equality constraints using a sequential quadratic programming (SQP) method. The minimization must be performed while observing only noisy evaluations of the objective and constraint functions. In order to obtain stability, the classical SQP method is modified by relaxing the standard Armijo line search based on the noise level in the functions, which is assumed to be known. Convergence theory is presented giving conditions under which the iterates converge to a neighborhood of the solution characterized by the noise level and the problem conditioning. The analysis assumes that the SQP algorithm does not require regularization or trust regions. Numerical experiments indicate that the relaxed line search improves the practical performance of the method on problems involving uniformly distributed noise, compared to a standard line search. [ABSTRACT FROM AUTHOR]

We develop a decomposition method based on the augmented Lagrangian framework to solve a broad family of semidefinite programming problems, possibly with nonlinear objective functions, nonsmooth regularization, and general linear equality/inequality constraints. In particular, the positive semidefinite variable along with a group of linear constraints can be transformed into a variable on a smooth manifold via matrix factorization. The nonsmooth regularization and other general linear constraints are handled by the augmented Lagrangian method. Therefore, each subproblem can be solved by a semismooth Newton method on a manifold. Theoretically, we show that the first and second-order necessary optimality conditions for the factorized subproblem are also sufficient for the original subproblem under certain conditions. Convergence analysis is established for the Riemannian subproblem and the augmented Lagrangian method. Extensive numerical experiments on large-scale semidefinite programming problems such as max-cut, nearest correlation estimation, clustering, and sparse principal component analysis demonstrate the strength of our proposed method compared to other state-of-the-art methods. [ABSTRACT FROM AUTHOR]

Point set matching problems can be handled by optimal transport. The mechanism behind it is that optimal transport recovers the point-to-point correspondence associated with the least curl deformation. Optimal transport is a special form of linear programming with dense constraints. Linear programming can be handled by interior point methods, provided that the involved ill-conditioned Hessians can be computed accurately. Matrix balancing has been employed to compute optimal transport under entropy regularization approaches. The solution quality relies on two factors: the accuracy of matrix balancing and the boundedness of the dual vector. High accurate matrix balancing is achieved by the application of Newton methods on a sequence of matrices along a central path. In this work, we apply sparse support constraints to matrix-balancing based interior point methods, in which the sparse set fulfilling total support is iteratively updated to truncate the domain of the transport plan. The total support condition is one crucial condition, which guarantees the existence of matrix balancing as well as the boundedness of the dual vector. [ABSTRACT FROM AUTHOR]