Error‐free transformation of matrix multiplication with a posteriori validation

In this study, we examine the accurate matrix multiplication in floating‐point arithmetic. We demonstrate the error‐free transformations of matrix multiplication using high performance basic linear algebra subprograms. These transformations can be applied to accurately compute the product of two matrices using floating‐point entries. A key technique for this calculation is error‐free splitting in floating‐point matrices. In this study, we improve upon our previous method by a posteriori validation using floating‐point exception. In the method, we utilize the presence of overflow in a positive manner for detecting whether rounding error occurs. If overflow occurs, the result contains some exceptional values such as and NaN , that is, the method fails by necessity. Otherwise, we can confirm that no rounding error occurs in the process. Therefore, reducing the possibility of overflow is important. The numerical results suggest that the proposed algorithm provides more accurate results compared with the original algorithm. Moreover, for the product of n × n matrices, when , the new algorithm reduces the computing time for error‐free transformation by an average of 20 % and up to 30 % compared with the original algorithm. Furthermore, the new algorithm can be used when matrix multiplication is performed using divide‐and‐conquer methods. Copyright © 2016 John Wiley & Sons, Ltd.     

A fast method to block-diagonalize a Hankel matrix

In this paper, we consider an approximate block diagonalization algorithm of an n×n real Hankel matrix in which the successive transformation matrices are upper triangular Toeplitz matrices, and propose a new fast approach to compute the factorization in O(n ²) operations. This method consists on using the revised Bini method (Lin et al., Theor Comp Sci 315: 511–523, 2004). To motivate our approach, we also propose an approximate factorization variant of the customary fast method based on Schur complementation adapted to the n×n real Hankel matrix. All algorithms have been implemented in Matlab and numerical results are included to illustrate the effectiveness of our approach.     

The bilinear complexity and practical algorithms for matrix multiplication

A method for deriving bilinear algorithms for matrix multiplication is proposed. New estimates for the bilinear complexity of a number of problems of the exact and approximate multiplication of rectangular matrices are obtained. In particular, the estimate for the boundary rank of multiplying 3 × 3 matrices is improved and a practical algorithm for the exact multiplication of square n × n matrices is proposed. The asymptotic arithmetic complexity of this algorithm is O(n ^2.7743).     

Fast matrix multiplication is stable

We perform forward error analysis for a large class of recursive matrix multiplication algorithms in the spirit of Bini and Lotti [Numer. Math. 36:63–72, 1980]. As a consequence of our analysis, we show that the exponent of matrix multiplication (the optimal running time) can be achieved by numerically stable algorithms. We also show that new group-theoretic algorithms proposed in Cohn and Umans [Foundations of Computer Science, 44th Annual IEEE Symposium, pp. 438–449, 2003] and Cohn et al. [Foundations of Computer Science, 46th Annual IEEE Symposium, pp. 379–388, 2005] are all included in the class of algorithms to which our analysis applies, and are therefore numerically stable. We perform detailed error analysis for three specific fast group-theoretic algorithms. J. Demmel acknowledges support of NSF under grants CCF-0444486, ACI-00090127, CNS-0325873 and of DOE under grant DE-FC02-01ER25478. I. Dumitriu acknowledges support of the Miller Institute for Basic Research in Science. R. Kleinberg is supported by an NSF Mathematical Sciences Postdoctoral Research Fellowship     

Fast algorithms for floating-point interval matrix multiplication

We discuss several methods for real interval matrix multiplication. First, earlier studies of fast algorithms for interval matrix multiplication are introduced: naive interval arithmetic, interval arithmetic by midpoint-radius form by Oishi-Rump and its fast variant by Ogita-Oishi. Next, three new and fast algorithms are developed. The proposed algorithms require one, two or three matrix products, respectively. The point is that our algorithms quickly predict which terms become dominant radii in interval computations. We propose a hybrid method to predict which algorithm is suitable for optimizing performance and width of the result. Numerical examples are presented to show the efficiency of the proposed algorithms.     

Estimates of the trace of the inverse of a symmetric matrix using the modified Chebyshev algorithm

In this paper we study how to compute an estimate of the trace of the inverse of a symmetric matrix by using Gauss quadrature and the modified Chebyshev algorithm. As auxiliary polynomials we use the shifted Chebyshev polynomials. Since this can be too costly in computer storage for large matrices we also propose to compute the modified moments with a stochastic approach due to Hutchinson (Commun Stat Simul 18:1059–1076, 1989). In memory of Gene H. Golub.     

Symmetric centrosymmetric matrix–vector multiplication

We present a method for the multiplication of an arbitrary vector by a symmetric centrosymmetric matrix, requiring floating-point operations, rather than the 2n² operations needed in the case of an arbitrary matrix. Combining this method with Trench's algorithm for Toeplitz matrix inversion yields a method for solving Toeplitz systems with the same complexity as Levinson's algorithm.     

A practical algorithm for faster matrix multiplication

The purpose of this paper is to present an algorithm for matrix multiplication based on a formula discovered by Pan [7]. For matrices of order up to 10 000, the nearly optimum tuning of the algorithm results in a rather clear non‐recursive one‐ or two‐level structure with the operation count comparable to that of the Strassen algorithm [9]. The algorithm takes less workspace and has better numerical stability as compared to the Strassen algorithm, especially in Winograd's modification [2]. Moreover, its clearer and more flexible structure is potentially more suitable for efficient implementation on modern supercomputers. Copyright © 1999 John Wiley & Sons, Ltd.     

Preconditioned conjugate gradient algorithms for nonconvex problems with box constraints

The paper describes new conjugate gradient algorithms for large scale nonconvex problems with box constraints. In order to speed up convergence the algorithms employ scaling matrices which transform the space of original variables into the space in which Hessian matrices of the problem’s functionals have more clustered eigenvalues. This is done by applying limited memory BFGS updating matrices. Once the scaling matrix is calculated, the next few conjugate gradient iterations are performed in the transformed space. The box constraints are treated efficiently by the projection. We also present a limited memory quasi-Newton method which is a special version of our general algorithm. The presented algorithms have strong global convergence properties, in particular they identify constraints active at a solution in a finite number of iterations. We believe that they are competitive to the L-BFGS-B method and present some numerical results which support our claim.     

Interpolation algorithm of Leverrier–Faddev type for polynomial matrices

We investigated an interpolation algorithm for computing outer inverses of a given polynomial matrix, based on the Leverrier–Faddeev method. This algorithm is a continuation of the finite algorithm for computing generalized inverses of a given polynomial matrix, introduced in [11]. Also, a method for estimating the degrees of polynomial matrices arising from the Leverrier–Faddeev algorithm is given as the improvement of the interpolation algorithm. Based on similar idea, we introduced methods for computing rank and index of polynomial matrix. All algorithms are implemented in the symbolic programming language MATHEMATICA , and tested on several different classes of test examples.     

Generalized matrix inversion is not harder than matrix multiplication

Starting from the Strassen method for rapid matrix multiplication and inversion as well as from the recursive Cholesky factorization algorithm, we introduced a completely block recursive algorithm for generalized Cholesky factorization of a given symmetric, positive semi-definite matrix A∈R^n×n$A\in {\mathbb{R}}^{n\times n}$. We used the Strassen method for matrix inversion together with the recursive generalized Cholesky factorization method, and established an algorithm for computing generalized {2,3}$\{2,3\}$ and {2,4}$\{2,4\}$ inverses. Introduced algorithms are not harder than the matrix–matrix multiplication.     

Deterministic computation of the characteristic polynomial in the time of matrix multiplication

This paper describes an algorithm which computes the characteristic polynomial of a matrix over a field within the same asymptotic complexity, up to constant factors, as the multiplication of two square matrices. Previously, this was only achieved by resorting to genericity assumptions or randomization techniques, while the best known complexity bound with a general deterministic algorithm was obtained by Keller-Gehrig in 1985 and involves logarithmic factors. Our algorithm computes more generally the determinant of a univariate polynomial matrix in reduced form, and relies on new subroutines for transforming shifted reduced matrices into shifted weak Popov matrices, and shifted weak Popov matrices into shifted Popov matrices.     

m-Fold Hypergeometric Solutions of Linear Recurrence Equations Revisited

We present two algorithms to compute m-fold hypergeometric solutions of linear recurrence equations for the classical shift case and for the q-case, respectively. The first is an m-fold generalization and q-generalization of the algorithm by van Hoeij (Appl Algebra Eng Commun Comput 17:83–115, 2005; J. Pure Appl Algebra 139:109–131, 1998) for recurrence equations. The second is a combination of an improved version of the algorithms by Petkovšek (Discrete Math 180:3–22, 1998; J Symb Comput 14(2–3):243–264, 1992) for recurrence and q-recurrence equations and the m-fold algorithm from Petkovšek and Salvy (ISSAC 1993 Proceedings, pp 27–33, 1993) for recurrence equations. We will refer to the classical algorithms as van Hoeij or Petkovšek respectively. To formulate our ideas, we first need to introduce an adapted version of an m-fold Newton polygon and its characteristic polynomials for the classical case and q-case, and to prove the important properties in this case. Using the data from the Newton polygon, we are able to present efficient m-fold versions of the van Hoeij and Petkovšek algorithms for the classical shift case and for the q-case, respectively. Furthermore, we show how one can use the Newton polygon and our characteristic polynomials to conclude for which m ? \mathbbN{m\in \mathbb{N}} there might be an m-fold hypergeometric solution at all. Again by using the information obtained from the Newton polygon, the presentation of the q-Petkovšek algorithm can be simplified and streamlined. Finally, we give timings for the ‘classical’ q-Petkovšek, our q-van Hoeij and our modified q-Petkovšek algorithm on some classes of problems and we present a Maple implementation of the m-fold algorithms for the q-case.     

Combinatorial Sublinear-Time Fourier Algorithms

We study the problem of estimating the best k term Fourier representation for a given frequency sparse signal (i.e., vector) A of length N≫k. More explicitly, we investigate how to deterministically identify k of the largest magnitude frequencies of [^(A)]\hat{\mathbf{A}} , and estimate their coefficients, in polynomial(k,log N) time. Randomized sublinear-time algorithms which have a small (controllable) probability of failure for each processed signal exist for solving this problem (Gilbert et al. in ACM STOC, pp. 152–161, 2002; Proceedings of SPIE Wavelets XI, 2005). In this paper we develop the first known deterministic sublinear-time sparse Fourier Transform algorithm which is guaranteed to produce accurate results. As an added bonus, a simple relaxation of our deterministic Fourier result leads to a new Monte Carlo Fourier algorithm with similar runtime/sampling bounds to the current best randomized Fourier method (Gilbert et al. in Proceedings of SPIE Wavelets XI, 2005). Finally, the Fourier algorithm we develop here implies a simpler optimized version of the deterministic compressed sensing method previously developed in (Iwen in Proc. of ACM-SIAM Symposium on Discrete Algorithms (SODA’08), 2008).     

Computing a matrix symmetrizer exactly using modified multiple modulus residue arithmetic

A symmetric solution X satisfying the matrix equation XA = A^tX is called a symmetrizer of the matrix A. A general algorithm to compute a matrix symmetrizer is obtained. A new multiple-modulus residue arithmetic called floating-point modular arithmetic is described and implemented on the algorithm to compute an error-free matrix symmetrizer.     

Jacobi spectral Galerkin method for elliptic Neumann problems

This paper is concerned with fast spectral-Galerkin Jacobi algorithms for solving one- and two-dimensional elliptic equations with homogeneous and nonhomogeneous Neumann boundary conditions. The paper extends the algorithms proposed by Shen (SIAM J Sci Comput 15:1489–1505, 1994) and Auteri et al. (J Comput Phys 185:427–444, 2003), based on Legendre polynomials, to Jacobi polynomials with arbitrary α and β. The key to the efficiency of our algorithms is to construct appropriate basis functions with zero slope at the endpoints, which lead to systems with sparse matrices for the discrete variational formulations. The direct solution algorithm developed for the homogeneous Neumann problem in two-dimensions relies upon a tensor product process. Nonhomogeneous Neumann data are accounted for by means of a lifting. Numerical results indicating the high accuracy and effectiveness of these algorithms are presented.     

An algorithm for automatically selecting a suitable verification method for linear systems

Several methods have been proposed to calculate a rigorous error bound of an approximate solution of a linear system by floating-point arithmetic. These methods are called ‘verification methods’. Applicable range of these methods are different. It depends mainly on the condition number and the dimension of the coefficient matrix whether such methods succeed to work or not. In general, however, the condition number is not known in advance. If the dimension or the condition number is large to some extent, then Oishi–Rump’s method, which is known as the fastest verification method for this purpose, may fail. There are more robust verification methods whose computational cost is larger than the Oishi–Rump’s one. It is not so efficient to apply such robust methods to well-conditioned problems. The aim of this paper is to choose a suitable verification method whose computational cost is minimum to succeed. First in this paper, four fast verification methods for linear systems are briefly reviewed. Next, a compromise method between Oishi–Rump’s and Ogita–Oishi’s one is developed. Then, an algorithm which automatically and efficiently chooses an appropriate verification method from five verification methods is proposed. The proposed algorithm does as much work as necessary to calculate error bounds of approximate solutions of linear systems. Finally, numerical results are presented.     

A Framework for the Construction of Self-replicating Tilings

Dekking (Adv. Math. 44:78–104, 1982; J. Comb. Theory Ser. A 32:315–320, 1982) provided an important method to compute the boundaries of lattice rep-tiles as a ‘recurrent set’ on a free group of a finite alphabet. That is, those tilings are generated by lattice translations of a single tile, and there is an expanding linear map that carries tiles to unions of tiles. The boundary of the tile is identified with a sequence of words in the alphabet obtained from an expanding endomorphism (substitution) on the alphabet. In this paper, Dekking’s construction is generalized to address tilings with more than one tile, and to have the elements of the tilings be generated by both translation and rotations. Examples that fall within the scope of our main result include self-replicating multi-tiles, self-replicating tiles for crystallographic tilings and aperiodic tilings.     

A hybrid method for quantum global optimization

This paper gives a quantum algorithm for global optimization. The heart of such approaches employ Grover’s database search (1996; Phys Rev Lett 79(23):4709–4712, 1997a; 79(2):325–328, 1997b). Chi and Kim (1998) show that when the phases of the generalized Grover database search operator are optimally chosen, it is capable of finding a solution by a single query. To apply this method to global optimization requires knowledge of the number of marked points m to calculate the optimal phases, but this value is seldom known. This paper focuses on overcoming this hurdle by showing that an estimate of the optimal phases can be found and used to replace the optimal phases while maintaining a high probability of finding a solution. Merging this finding with a recently discovered dynamic quantum global optimization algorithm (BBW2D) that reduces the problem to finding successively improving regions using Grover’s search, we present a hybrid method that improves the efficiency and reduces the variance of the search algorithm when empirically compared to other existing quantum search algorithms.