Search

The success of large language models (LLMs) for time series has been demonstrated in previous work. Utilizing a symbolic time series representation, one can efficiently bridge the gap between LLMs and time series. However, the remaining challenge is to exploit the semantic information hidden in time series by using symbols or existing tokens of LLMs, while aligning the embedding space of LLMs according to the hidden information of time series. The symbolic time series approximation (STSA) method called adaptive Brownian bridge-based symbolic aggregation (ABBA) shows outstanding efficacy in preserving salient time series features by modeling time series patterns in terms of amplitude and period while using existing tokens of LLMs. In this paper, we introduce a method, called LLM-ABBA, that integrates ABBA into large language models for various downstream time series tasks. By symbolizing time series, LLM-ABBA compares favorably to the recent state-of-the-art (SOTA) in UCR and three medical time series classification tasks. Meanwhile, a fixed-polygonal chain trick in ABBA is introduced to \kc{avoid obvious drifting} during prediction tasks by significantly mitigating the effects of cumulative error arising from misused symbols during the transition from symbols to numerical values. In time series regression tasks, LLM-ABBA achieves the new SOTA on Time Series Extrinsic Regression (TSER) benchmarks. LLM-ABBA also shows competitive prediction capability compared to recent SOTA time series prediction results. We believe this framework can also seamlessly extend to other time series tasks.

Time series are ubiquitous in numerous science and engineering domains, e.g., signal processing, bioinformatics, and astronomy. Previous work has verified the efficacy of symbolic time series representation in a variety of engineering applications due to its storage efficiency and numerosity reduction. The most recent symbolic aggregate approximation technique, ABBA, has been shown to preserve essential shape information of time series and improve downstream applications, e.g., neural network inference regarding prediction and anomaly detection in time series. Motivated by the emergence of high-performance hardware which enables efficient computation for low bit-width representations, we present a new quantization-based ABBA symbolic approximation technique, QABBA, which exhibits improved storage efficiency while retaining the original speed and accuracy of symbolic reconstruction. We prove an upper bound for the error arising from quantization and discuss how the number of bits should be chosen to balance this with other errors. An application of QABBA with large language models (LLMs) for time series regression is also presented, and its utility is investigated. By representing the symbolic chain of patterns on time series, QABBA not only avoids the training of embedding from scratch, but also achieves a new state-of-the-art on Monash regression dataset. The symbolic approximation to the time series offers a more efficient way to fine-tune LLMs on the time series regression task which contains various application domains. We further present a set of extensive experiments performed across various well-established datasets to demonstrate the advantages of the QABBA method for symbolic approximation.

The block classical Gram--Schmidt (BCGS) algorithm and its reorthogonalized variant are widely-used methods for computing the economic QR factorization of block columns $X$ due to their lower communication cost compared to other approaches such as modified Gram--Schmidt and Householder QR. To further reduce communication, i.e., synchronization, there has been a long ongoing search for a variant of reorthogonalized BCGS variant that achieves $O(u)$ loss of orthogonality while requiring only \emph{one} synchronization point per block column, where $u$ represents the unit roundoff. Utilizing Pythagorean inner products and delayed normalization techniques, we propose the first provably stable one-synchronization reorthogonalized BCGS variant, demonstrating that it has $O(u)$ loss of orthogonality under the condition $O(u) \kappa^2(X) \leq 1/2$, where $\kappa(\cdot)$ represents the condition number. By incorporating one additional synchronization point, we develop a two-synchronization reorthogonalized BCGS variant which maintains $O(u)$ loss of orthogonality under the improved condition $O(u) \kappa(X) \leq 1/2$. An adaptive strategy is then proposed to combine these two variants, ensuring $O(u)$ loss of orthogonality while using as few synchronization points as possible under the less restrictive condition $O(u) \kappa(X) \leq 1/2$. As an example of where this adaptive approach is beneficial, we show that using the adaptive orthogonalization variant, $s$-step GMRES achieves a backward error comparable to $s$-step GMRES with BCGSI+, also known as BCGS2, both theoretically and numerically, but requires fewer synchronization points., Comment: 25 pages, 8. figures

Sketching-based preconditioners have been shown to accelerate the solution of dense least-squares problems with coefficient matrices having substantially more rows than columns. The cost of generating these preconditioners can be reduced by employing low precision floating-point formats for all or part of the computations. We perform finite precision analysis of a mixed precision algorithm that computes the $R$-factor of a QR factorization of the sketched coefficient matrix. Two precisions can be chosen and the analysis allows understanding how to set these precisions to exploit the potential benefits of low precision formats and still guarantee an effective preconditioner. If the nature of the least-squares problem requires a solution with a small forward error, then mixed precision iterative refinement (IR) may be needed. For ill-conditioned problems the GMRES-based IR approach can be used, but good preconditioner is crucial to ensure convergence. We theoretically show when the sketching-based preconditioner can guarantee that the GMRES-based IR reduces the relative forward error of the least-squares solution and the residual to the level of the working precision unit roundoff. Small numerical examples illustrate the analysis.

As computational machines are becoming larger and more complex, the probability of hardware failure rises. ``Silent errors'', or, bit flips, may not be immediately apparent but can cause detrimental effects to algorithm behavior. In this work, we examine an algorithm-based approach to silent error detection in the context of pipelined Krylov subspace methods, in particular, Pipe-PR-CG, for the solution of linear systems. Our approach is based on using finite precision error analysis to bound the differences between quantities which should be equal in exact arithmetic. Through inexpensive monitoring during the iteration, we can detect when these bounds are violated, which indicates that a silent error has occurred. We use this approach to develop a fault-tolerance variant and also suggest a strategy for dynamically adapting the detection criteria. Our numerical experiments demonstrate the effectiveness of our approach., Comment: 42 pages

Communication, i.e., data movement, is a critical bottleneck for the performance of classical Krylov subspace method solvers on modern computer architectures. Variants of these methods which avoid communication have been introduced, which, while equivalent in exact arithmetic, can be unstable in finite precision. In this work, we address the backward stability of s-step GMRES, also known as communication-avoiding GMRES. We present a framework for simplifying the analysis of s-step GMRES, which includes standard GMRES (s=1) as a special case, by isolating the effects of rounding errors in the QR factorization and the solution of the least squares problem. Using this framework, we analyze s-step GMRES with popular block orthogonalization methods: block modified Gram--Schmidt and reorthogonalized block classical Gram--Schmidt algorithms. An example illustrates the resulting instability of s-step GMRES when paired with the classical s-step Arnoldi process and shows the limitations of popular strategies for resolving this instability. To address this issue, we propose a modified Arnoldi process that allows for much larger block size s while maintaining satisfactory accuracy, as confirmed by our numerical experiments., Comment: 32 pages, 9 figures

Interest in communication-avoiding orthogonalization schemes for high-performance computing has been growing recently. This manuscript addresses open questions about the numerical stability of various block classical Gram-Schmidt variants that have been proposed in the past few years. An abstract framework is employed, the flexibility of which allows for new rigorous bounds on the loss of orthogonality in these variants. We first analyze a generalization of (reorthogonalized) block classical Gram-Schmidt and show that a "strong" intrablock orthogonalization routine is only needed for the very first block in order to maintain orthogonality on the level of the unit roundoff. Then, using this variant, which has four synchronization points per block column, we remove the synchronization points one at a time and analyze how each alteration affects the stability of the resulting method. Our analysis shows that the variant requiring only one synchronization per block column cannot be guaranteed to be stable in practice, as stability begins to degrade with the first reduction of synchronization points. Our analysis of block methods also provides new theoretical results for the single-column case. In particular, it is proven that DCGS2 from [Bielich, D. et al. Par. Comput. 112 (2022)] and CGS-2 from [\'{S}wirydowicz, K. et al, Num. Lin. Alg. Appl. 28 (2021)] are as stable as Householder QR. Numerical examples from the BlockStab toolbox are included throughout, to help compare variants and illustrate the effects of different choices of intraorthogonalization subroutines.

Hierarchical matrix computations have attracted significant attention in the science and engineering community as exploiting data-sparse structures can significantly reduce the computational complexity of many important kernels. One particularly popular option within this class is the Hierarchical Off-Diagonal Low-Rank (HODLR) format. In this paper, we show that the off-diagonal blocks of HODLR matrices that are approximated by low-rank matrices can be represented in low precision without degenerating the quality of the overall approximation (with the error growth bounded by a factor of $2$). We also present an adaptive-precision scheme for constructing and storing HODLR matrices, and we prove that the use of mixed precision does not compromise the numerical stability of the resulting HODLR matrix--vector product and LU factorization. That is, the resulting error in these computations is not significantly greater than the case where we use one precision (say, double) for constructing and storing the HODLR matrix. Our analyses further give insight on how one must choose the working precision in HODLR matrix computations relative to the approximation error in order to not observe the effects of finite precision. Intuitively, when a HODLR matrix is subject to a high degree of approximation error, subsequent computations can be performed in a lower precision without detriment. We demonstrate the validity of our theoretical results through a range of numerical experiments.

The $k$-means algorithm is one of the most popular and critical techniques in data mining and machine learning, and it has achieved significant success in numerous science and engineering domains. Computing $k$-means to a global optimum is NP-hard in Euclidean space, yet there are a variety of efficient heuristic algorithms, such as Lloyd's algorithm, that converge to a local optimum with superpolynomial complexity in the worst case. Motivated by the emergence and prominence of mixed precision capabilities in hardware, a current trend is to develop low and mixed precision variants of algorithms in order to improve the runtime and energy consumption. In this paper we study the numerical stability of Lloyd's $k$-means algorithm, and, in particular, we confirm the stability of the widely used distance computation formula. We propose a mixed-precision framework for $k$-means computation and investigate the effects of low-precision distance computation within the framework. Through extensive simulations on various data clustering and image segmentation tasks, we verify the applicability and robustness of the mixed precision $k$-means method. We find that, in $k$-means computation, normalized data is more tolerant to the reduction of precision in the distance computation, while for nonnormalized data more care is needed in the use of reduced precision, mainly to avoid overflow. Our study demonstrates the potential for the use of mixed precision to accelerate the $k$-means computation and offers some insights into other distance-based machine learning methods.

Various approaches to iterative refinement (IR) for least-squares problems have been proposed in the literature and it may not be clear which approach is suitable for a given problem. We consider three approaches to IR for least-squares problems when two precisions are used and review their theoretical guarantees, known shortcomings and when the method can be expected to recognize that the correct solution has been found, and extend uniform precision analysis for an IR approach based on the semi-normal equations to the two-precision case. We focus on the situation where it is desired to refine the solution to the working precision level. It is shown that the IR methods exhibit different sensitivities to the conditioning of the problem and the size of the least-squares residual, which should be taken into account when choosing the IR approach. We also discuss a new approach that is based on solving multiple least-squares problems.

Block classical Gram-Schmidt (BCGS) is commonly used for orthogonalizing a set of vectors $X$ in distributed computing environments due to its favorable communication properties relative to other orthogonalization approaches, such as modified Gram-Schmidt or Householder. However, it is known that BCGS (as well as recently developed low-synchronization variants of BCGS) can suffer from a significant loss of orthogonality in finite-precision arithmetic, which can contribute to instability and inaccurate solutions in downstream applications such as $s$-step Krylov subspace methods. A common solution to improve the orthogonality among the vectors is reorthogonalization. Focusing on the "Pythagorean" variant of BCGS, introduced in [E. Carson, K. Lund, & M. Rozlo\v{z}n\'{i}k. SIAM J. Matrix Anal. Appl. 42(3), pp. 1365--1380, 2021], which guarantees an $O(\varepsilon)\kappa^2(X)$ bound on the loss of orthogonality as long as $O(\varepsilon)\kappa^2(X)<1$, where $\varepsilon$ denotes the unit roundoff, we introduce and analyze two reorthogonalized Pythagorean BCGS variants. These variants feature favorable communication properties, with asymptotically two synchronization points per block column, as well as an improved $O(\varepsilon)$ bound on the loss of orthogonality. Our bounds are derived in a general fashion to additionally allow for the analysis of mixed-precision variants. We verify our theoretical results with a panel of test matrices and experiments from a new version of the \texttt{BlockStab} toolbox.

With the recent emergence of mixed precision hardware, there has been a renewed interest in its use for solving numerical linear algebra problems fast and accurately. The solution of least squares (LS) problems $\min_x\|b-Ax\|_2$, where $A \in \mathbb{R}^{m\times n}$, arise in numerous application areas. Overdetermined standard least squares problems can be solved by using mixed precision within the iterative refinement method of Bj\"{o}rck, which transforms the least squares problem into an $(m+n)\times(m+n)$ ''augmented'' system. It has recently been shown that mixed precision GMRES-based iterative refinement can also be used, in an approach termed GMRES-LSIR. In practice, we often encounter types of least squares problems beyond standard least squares, including weighted least squares (WLS), $\min_x\|D^{1/2}(b-Ax)\|_2$, where $D^{1/2}$ is a diagonal matrix of weights. In this paper, we discuss a mixed precision FGMRES-WLSIR algorithm for solving WLS problems using two different preconditioners., Comment: 12 pages

Hardware trends have motivated the development of mixed precision algo-rithms in numerical linear algebra, which aim to decrease runtime while maintaining acceptable accuracy. One recent development is the development of an adaptive precision sparse matrix-vector produce routine, which may be used to accelerate the solution of sparse linear systems by iterative methods. This approach is also applicable to the application of inexact preconditioners, such as sparse approximate inverse preconditioners used in Krylov subspace methods. In this work, we develop an adaptive precision sparse approximate inverse preconditioner and demonstrate its use within a five-precision GMRES-based iterative refinement method. We call this algorithm variant BSPAI-GMRES-IR. We then analyze the conditions for the convergence of BSPAI-GMRES-IR, and determine settings under which BSPAI-GMRES-IR will produce similar backward and forward errors as the existing SPAI-GMRES-IR method, the latter of which does not use adaptive precision in preconditioning. Our numerical experiments show that this approach can potentially lead to a reduction in the cost of storing and applying sparse approximate inverse preconditioners, although a significant reduction in cost may comes at the expense of increasing the number of GMRES iterations required for convergence.

The multigrid V-cycle method is a popular method for solving systems of linear equations. It computes an approximate solution by using smoothing on fine levels and solving a system of linear equations on the coarsest level. Solving on the coarsest level depends on the size and difficulty of the problem. If the size permits, it is typical to use a direct method based on LU or Cholesky decomposition. In settings with large coarsest-level problems, approximate solvers such as iterative Krylov subspace methods, or direct methods based on low-rank approximation, are often used. The accuracy of the coarsest-level solver is typically determined based on the experience of the users with the concrete problems and methods. In this paper we present an approach to analyzing the effects of approximate coarsest-level solves on the convergence of the V-cycle method for symmetric positive definite problems. Using these results, we derive coarsest-level stopping criterion through which we may control the difference between the approximation computed by a V-cycle method with approximate coarsest-level solver and the approximation which would be computed if the coarsest-level problems were solved exactly. The coarsest-level stopping criterion may thus be set up such that the V-cycle method converges to a chosen finest-level accuracy in (nearly) the same number of V-cycle iterations as the V-cycle method with exact coarsest-level solver. We also utilize the theoretical results to discuss how the convergence of the V-cycle method may be affected by the choice of a tolerance in a coarsest-level stopping criterion based on the relative residual norm., Comment: 26 pages, 8 figures

With the recent emergence of mixed precision hardware, there has been a renewed interest in its use for solving numerical linear algebra problems fast and accurately. The solution of total least squares problems, i.e., solving $\min_{E,r} \| [E, r]\|_F$ subject to $(A+E)x=b+r$, arises in numerous applications. Solving this problem requires finding the smallest singular value and corresponding right singular vector of $[A,b]$, which is challenging when $A$ is large and sparse. An efficient algorithm for this case due to Bj\"{o}rck et al. [SIAM J. Matrix Anal. Appl. 22(2), 2000], called RQI-PCGTLS, is based on Rayleigh quotient iteration coupled with the preconditioned conjugate gradient method. We develop a mixed precision variant of this algorithm, RQI-PCGTLS-MP, in which up to three different precisions can be used. We assume that the lowest precision is used in the computation of the preconditioner, and give theoretical constraints on how this precision must be chosen to ensure stability. In contrast to standard least squares, for total least squares, the resulting constraint depends not only on the matrix $A$, but also on the right-hand side $b$. We perform a number of numerical experiments on model total least squares problems used in the literature, which demonstrate that our algorithm can attain the same accuracy as RQI-PCGTLS albeit with a potential convergence delay due to the use of low precision. Performance modeling shows that the mixed precision approach can achieve up to a $4\times$ speedup depending on the size of the matrix and the number of Rayleigh quotient iterations performed., Comment: 20 pages

We consider the split-preconditioned FGMRES method in a mixed precision framework, in which four potentially different precisions can be used for computations with the coefficient matrix, application of the left preconditioner, application of the right preconditioner, and the working precision. Our analysis is applicable to general preconditioners. We obtain bounds on the backward and forward errors in split-preconditioned FGMRES. Our analysis further provides insight into how the various precisions should be chosen; under certain assumptions, a suitable selection guarantees a backward error on the order of the working precision.

When the CG method for solving linear algebraic systems was formulated about 70 years ago by Lanczos, Hestenes, and Stiefel, it was considered an iterative process possessing a mathematical finite termination property. CG was placed into a rich mathematical context, including links with Gauss quadrature and continued fractions. The optimality property of CG was described via a normalized weighted polynomial least squares approximation to zero. This highly nonlinear problem explains the adaptation of CG iterates to the given data. Karush and Hayes immediately considered CG in infinite dimensional Hilbert spaces and investigated its superlinear convergence. Since then, the view of CG and other Krylov subspace methods has changed. Today these methods are primarily used as computational tools, and their behavior is typically characterized using linear upper bounds or heuristics based on clustering of eigenvalues. Such simplifications limit the mathematical understanding and also negatively affect their practical application. This paper offers a different perspective. Focusing on CG and GMRES, it presents mathematically important and practically relevant phenomena that uncover their behavior through a discussion of computed examples. These examples provide an easily accessible approach that enables understanding of the methods, while pointers to more detailed analyses in the literature are given. This approach allows readers to choose the level of depth and thoroughness appropriate for their intentions. Some of the points made in this paper illustrate well known facts. Others challenge mainstream views and explain existing misunderstandings. Several points refer to recent results leading to open problems. We consider CG and GMRES crucially important for the mathematical understanding, further development, and practical applications also of other Krylov subspace methods., Comment: 56 pages

Using lower precision in algorithms can be beneficial in terms of reducing both computation and communication costs. Motivated by this, we aim to further the state-of-the-art in developing and analyzing mixed precision variants of iterative methods. In this work, we focus on the block variant of low-synchronization classical Gram-Schmidt with reorthogonalization, which we call BCGSI+LS. We demonstrate that the loss of orthogonality produced by this orthogonalization scheme can exceed $O(u)\kappa(\mathcal{X})$, where $u$ is the unit roundoff and $\kappa(\mathcal{X})$ is the condition number of the matrix to be orthogonalized, and thus we can not in general expect this to result in a backward stable block GMRES implementation. We then develop a mixed precision variant of this algorithm, called BCGSI+LS-MP, which uses higher precision in certain parts of the computation. We demonstrate experimentally that for a number of challenging test problems, our mixed precision variant successfully maintains a loss of orthogonality below $O(u)\kappa(\mathcal{X})$. This indicates that we can achieve a backward stable block GMRES algorithm that requires only one synchronization per iteration., Comment: 9 pages

Low rank matrix approximations appear in a number of scientific computing applications. We consider the Nystr\"{o}m method for approximating a positive semidefinite matrix $A$. In the case that $A$ is very large or its entries can only be accessed once, a single-pass version may be necessary. In this work, we perform a complete rounding error analysis of the single-pass Nystr\"{o}m method in two precisions, where the computation of the expensive matrix product with $A$ is assumed to be performed in the lower of the two precisions. Our analysis gives insight into how the sketching matrix and shift should be chosen to ensure stability, implementation aspects which have been commented on in the literature but not yet rigorously justified. We further develop a heuristic to determine how to pick the lower precision, which confirms the general intuition that the lower the desired rank of the approximation, the lower the precision we can use without detriment. We also demonstrate that our mixed precision Nystr\"{o}m method can be used to inexpensively construct limited memory preconditioners for the conjugate gradient method and derive a bound the condition number of the resulting preconditioned coefficient matrix. We present numerical experiments on a set of matrices with various spectral decays and demonstrate the utility of our mixed precision approach.

The GMRES algorithm of Saad and Schultz (1986) is an iterative method for approximately solving linear systems $A{\bf x}={\bf b}$, with initial guess ${\bf x}_0$ and residual ${\bf r}_0 = {\bf b} - A{\bf x}_0$. The algorithm employs the Arnoldi process to generate the Krylov basis vectors (the columns of $V_k$). It is well known that this process can be viewed as a $QR$ factorization of the matrix $B_k = [\: {\bf r}_0, AV_k\:]$ at each iteration. Despite an ${O}(\epsilon)\kappa(B_k)$ loss of orthogonality, for unit roundoff $\epsilon$ and condition number $\kappa$, the modified Gram-Schmidt formulation was shown to be backward stable in the seminal paper by Paige et al. (2006). We present an iterated Gauss-Seidel formulation of the GMRES algorithm (IGS-GMRES) based on the ideas of Ruhe (1983) and \'{S}wirydowicz et al. (2020). IGS-GMRES maintains orthogonality to the level ${O}(\epsilon)\kappa(B_k)$ or ${O}(\epsilon)$, depending on the choice of one or two iterations; for two Gauss-Seidel iterations, the computed Krylov basis vectors remain orthogonal to working precision and the smallest singular value of $V_k$ 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 non-normal systems., Comment: Updates to funding

With the commercial availability of mixed precision hardware, mixed precision GMRES-based iterative refinement schemes have emerged as popular approaches for solving sparse linear systems. Existing analyses of these approaches, however, are based on using full LU factorizations to construct preconditioners for use within GMRES in each refinement step. In practical applications, inexact preconditioning techniques, such as incomplete LU or sparse approximate inverses, are often used for performance reasons. In this work, we investigate the use of sparse approximate inverse preconditioners based on Frobenius norm minimization within GMRES-based iterative refinement. We analyze the computation of sparse approximate inverses in finite precision and derive constraints under which user-specified stopping criteria will be satisfied. We then analyze the behavior of and convergence constraints for a five-precision GMRES-based iterative refinement scheme that uses sparse approximate inverse preconditioning, which we call SPAI-GMRES-IR. Our numerical experiments confirm the theoretical analysis and illustrate the resulting tradeoffs between preconditioner sparsity and GMRES-IR convergence rate., Comment: 23 pages

With the emergence of mixed precision capabilities in hardware, iterative refinement schemes for solving linear systems $Ax=b$ have recently been revisited and reanalyzed in the context of three or more precisions. These new analyses show that under certain constraints on condition number, the LU factorization of the matrix can be computed in low precision without affecting the final accuracy. Another promising technique is GMRES-based iterative refinement, which, in contrast to the standard approach, use GMRES preconditioned by the low-precision triangular factors to solve for the approximate solution update in each refinement step. This more accurate solution method extends the range of problems which can be solved with a given combination of precisions. However, in certain settings, GMRES may require too many iterations per refinement step, making it potentially more expensive than simply recomputing the LU factors in a higher precision. Krylov subspace recycling is a well-known technique for reusing information across sequential invocations of a Krylov subspace method on systems with the same or a slowly changing coefficient matrix. In this work, we incorporate the idea of Krylov subspace recycling into a mixed precision GMRES-based iterative refinement solver. The insight is that in each refinement step, we call preconditioned GMRES on a linear system with the same coefficient matrix $A$, with only the right-hand side changing. In this way, the GMRES solves in subsequent refinement steps can be accelerated by recycling information obtained from the first step. We perform extensive numerical experiments on various random dense problems, Toeplitz problems (prolate matrices), and problems from real applications, which confirm the benefits of the recycling approach., Comment: 18 pages

Low precision arithmetic, in particular half precision floating point arithmetic, is now available in commercial hardware. Using lower precision can offer significant savings in computation and communication costs with proportional savings in energy. Motivated by this, there has been a renewed interest in mixed precision iterative refinement for solving linear systems $Ax=b$, and new variants of GMRES-based iterative refinement have been developed. Each particular variant with a given combination of precisions leads to different condition number-based constraints for convergence of the backward and forward errors, and each has different performance costs. The constraints for convergence given in the literature are, as an artifact of the analyses, often overly strict in practice, and thus could lead a user to select a more expensive variant when a less expensive one would have sufficed. In this work, we develop a multistage mixed precision iterative refinement solver which aims to combine existing mixed precision approaches to balance performance and accuracy and improve usability. For an initial combination of precisions, the algorithm begins with the least expensive approach and convergence is monitored via inexpensive computations with quantities produced during the iteration. If slow convergence or divergence is detected using particular stopping criteria, the algorithm switches to use a more expensive, but more reliable variant. A novel aspect of our approach is that, unlike existing implementations, our algorithm first attempts to use ``stronger'' solvers for the solution update before resorting to increasing the precision(s). In some scenarios, this can avoid the need to refactorize the matrix in higher precision. We perform extensive numerical experiments on random dense problems and problems from real applications which confirm the benefits of the multistage approach., Comment: 30 pages

Compared to the classical Lanczos algorithm, the $s$-step Lanczos variant has the potential to improve performance by asymptotically decreasing the synchronization cost per iteration. However, this comes at a cost. Despite being mathematically equivalent, the $s$-step variant is known to behave quite differently in finite precision, with potential for greater loss of accuracy and a decrease in the convergence rate relative to the classical algorithm. It has previously been shown that the errors that occur in the $s$-step version follow the same structure as the errors in the classical algorithm, but with the addition of an amplification factor that depends on the square of the condition number of the $O(s)-$dimensional Krylov bases computed in each outer loop. As the condition number of these $s$-step bases grows (in some cases very quickly) with $s$, this limits the parameter $s$ that can be chosen and thus limits the performance that can be achieved. In this work we show that if a select few computations in $s$-step Lanczos are performed in double the working precision, the error terms then depend only linearly on the conditioning of the $s$-step bases. This has the potential for drastically improving the numerical behavior of the algorithm with little impact on per-iteration performance. Our numerical experiments demonstrate the improved numerical behavior possible with the mixed precision approach, and also show that this improved behavior extends to the $s$-step CG algorithm in mixed precision., Comment: 34 pages

Block Gram-Schmidt algorithms serve as essential kernels in many scientific computing applications, but for many commonly used variants, a rigorous treatment of their stability properties remains open. This work provides a comprehensive categorization of block Gram-Schmidt algorithms, particularly those used in Krylov subspace methods to build orthonormal bases one block vector at a time. Known stability results are assembled, and new results are summarized or conjectured for important communication-reducing variants. Additionally, new block versions of low-synchronization variants are derived, and their efficacy and stability are demonstrated for a wide range of challenging examples. Numerical examples are computed with a versatile MATLAB package hosted at https://github.com/katlund/BlockStab, and scripts for reproducing all results in the paper are provided. Block Gram-Schmidt implementations in popular software packages are discussed, along with a number of open problems. An appendix containing all algorithms type-set in a uniform fashion is provided., Comment: 53 pages, 5 tables, 16 figures, 22 algorithms

Within the past years, hardware vendors have started designing low precision special function units in response to the demand of the Machine Learning community and their demand for high compute power in low precision formats. Also the server-line products are increasingly featuring low-precision special function units, such as the NVIDIA tensor cores in ORNL's Summit supercomputer providing more than an order of magnitude higher performance than what is available in IEEE double precision. At the same time, the gap between the compute power on the one hand and the memory bandwidth on the other hand keeps increasing, making data access and communication prohibitively expensive compared to arithmetic operations. To start the multiprecision focus effort, we survey the numerical linear algebra community and summarize all existing multiprecision knowledge, expertise, and software capabilities in this landscape analysis report. We also include current efforts and preliminary results that may not yet be considered "mature technology," but have the potential to grow into production quality within the multiprecision focus effort. As we expect the reader to be familiar with the basics of numerical linear algebra, we refrain from providing a detailed background on the algorithms themselves but focus on how mixed- and multiprecision technology can help improving the performance of these methods and present highlights of application significantly outperforming the traditional fixed precision methods., Comment: Technical report as a part of the Exascale computing project (ECP)

The adaptive $s$-step CG algorithm is a solver for sparse, symmetric positive definite linear systems designed to reduce the synchronization cost per iteration while still achieving a user-specified accuracy requirement. In this work, we improve the adaptive $s$-step conjugate gradient algorithm by use of iteratively updated estimates of the largest and smallest Ritz values, which give approximations of the largest and smallest eigenvalues of $A$, using a technique due to Meurant and Tich{\' y} [G. Meurant and P. Tich{\' y}, Numer. Algs. (2018), pp.~1--32]. The Ritz value estimates are used to dynamically update parameters for constructing Newton or Chebyshev polynomials so that the conditioning of the $s$-step bases can be continuously improved throughout the iterations. These estimates are also used to automatically set a variable related to the ratio of the sizes of the error and residual, which was previously treated as an input parameter. We show through numerical experiments that in many cases the new algorithm improves upon the previous adaptive $s$-step approach both in terms of numerical behavior and reduction in number of synchronizations., Comment: 25 pages, 4 figures

The standard implementation of the conjugate gradient algorithm suffers from communication bottlenecks on parallel architectures, due primarily to the two global reductions required every iteration. In this paper, we study conjugate gradient variants which decrease the runtime per iteration by overlapping global synchronizations, and in the case of pipelined variants, matrix-vector products. Through the use of a predict-and-recompute scheme, whereby recursively-updated quantities are first used as a predictor for their true values and then recomputed exactly at a later point in the iteration, these variants are observed to have convergence behavior nearly as good as the standard conjugate gradient implementation on a variety of test problems. We provide a rounding error analysis which provides insight into this observation. It is also verified experimentally that the variants studied do indeed reduce the runtime per iteration in practice and that they scale similarly to previously-studied communication-hiding variants. Finally, because these variants achieve good convergence without the use of any additional input parameters, they have the potential to be used in place of the standard conjugate gradient implementation in a range of applications., Comment: This material is based upon work supported by the NSF GRFP. Code for reproducing all figures and tables can be found here: https://github.com/tchen01/new_cg_variants

On modern large-scale parallel computers, the performance of Krylov subspace iterative methods is limited by global synchronization. This has inspired the development of $s$-step Krylov subspace method variants, in which iterations are computed in blocks of $s$, which can reduce the number of global synchronizations per iteration by a factor of $O(s)$. Although the $s$-step variants are mathematically equivalent to their classical counterparts, they can behave quite differently in finite precision depending on the parameter $s$. If $s$ is chosen too large, the $s$-step method can suffer a convergence delay and a decrease in attainable accuracy relative to the classical method. This makes it difficult for a potential user of such methods - the $s$ value that minimizes the time per iteration may not be the best $s$ for minimizing the overall time-to-solution, and further may cause an unacceptable decrease in accuracy. Towards improving the reliability and usability of $s$-step Krylov subspace methods, in this work we derive the \emph{adaptive $s$-step CG method}, a variable $s$-step CG method where in block $k$, the parameter $s_k$ is determined automatically such that a user-specified accuracy is attainable. The method for determining $s_k$ is based on a bound on growth of the residual gap within block $k$, from which we derive a constraint on the condition numbers of the computed $O(s_k)$-dimensional Krylov subspace bases. The computations required for determining the block size $s_k$ can be performed without increasing the number of global synchronizations per block. Our numerical experiments demonstrate that the adaptive $s$-step CG method is able to attain up to the same accuracy as classical CG while still significantly reducing the total number of global synchronizations.