Fast rectangular matrix multiplication and some applications

We study asymptotically fast multiplication algorithms for matrix pairs of arbitrary di- mensions, and optimize the exponents of their arithmetic complexity bounds. For a large class of input matrix pairs, we improve the known exponents. We also show some applications of our results:(i) we decrease from O(n~2 n~(1 o)(1)logq)to O(n~(1.9998) n~(1 o(1))logq)the known arithmetic complexity bound for the univariate polynomial factorization of degree n over a finite field with q elements; (ii) we decrease from 2.837 to 2.7945 the known exponent of the work and arithmetic processor bounds for fast deterministic(NC)parallel evaluation of the determinant, the characteristic polynomial, and the inverse of an n×n matrix, as well as for the solution to a nonsingular linear system of n equations; (iii)we decrease from O(m~(1.575)n)to O(m~(1.5356)n)the known bound for computing basic solutions to a linear programming problem with m constraints and n variables.     

Computing matrix functions solving coupled differential models

In this paper a modification of the method proposed in [E. Defez, L. Jódar, Some applications of Hermite matrix polynomials series expansions, Journal of Computational and Applied Mathematics 99 (1998) 105–117] for computing matrix sine and cosine based on Hermite matrix polynomial expansions is presented. An algorithm and illustrative examples demonstrate the performance of the new proposed method.     

Polynomial numerical hulls of matrix polynomials,II

In this article, we study some algebraic and geometrical properties of polynomial numerical hulls of matrix polynomials and joint polynomial numerical hulls of a finite family of matrices (possibly the coefficients of a matrix polynomial). Also, we study polynomial numerical hulls of basic A-factor block circulant matrices. These are block companion matrices of particular simple monic matrix polynomials. By studying the polynomial numerical hulls of the Kronecker product of two matrices, we characterize the polynomial numerical hulls of unitary basic A-factor block circulant matrices.     

Differentiation of generalized inverses for rational and polynomial matrices

Our basic motivation is a direct method for computing the gradient of the pseudo-inverse of well-conditioned system with respect to a scalar, proposed in [13] by Layton. In the present paper we combine the Layton’s method together with the representation of the Moore-Penrose inverse of one-variable polynomial matrix from [24] and developed an algorithm for computing the gradient of the Moore-Penrose inverse for one-variable polynomial matrix. Moreover, using the representation of various types of pseudo-inverses from [26], based on the Grevile’s partitioning method, we derive more general algorithms for computing {1}, {1, 3} and {1, 4} inverses of one-variable rational and polynomial matrices. Introduced algorithms are implemented in the programming language MATHEMATICA. Illustrative examples on analytical matrices are presented.     

Roundoff error analysis of fast DCT algorithms in fixed point arithmetic

Discrete cosine transforms (DCT) are essential tools in numerical analysis and digital signal processing. Processors in digital signal processing often use fixed point arithmetic. In this paper, we consider the numerical stability of fast DCT algorithms in fixed point arithmetic. The fast DCT algorithms are based on known factorizations of the corresponding cosine matrices into products of sparse, orthogonal matrices of simple structure. These algorithms are completely recursive, are easy to implement and use only permutations, scaling, butterfly operations, and plane rotations/rotation-reflections. In comparison with other fast DCT algorithms, these algorithms have low arithmetic costs. Using von Neumann–Goldstine’s model of fixed point arithmetic, we present a detailed roundoff error analysis for fast DCT algorithms in fixed point arithmetic. Numerical tests demonstrate the performance of our results.      

A note on the computation of the CP-rank

The purpose of this note is to address the computational question of determining whether or not a square nonnegative matrix (over the reals) is completely positive and finding its CP-rank when it is. We show that these questions can be resolved by finite algorithms and we provide (non-polynomial) complexity bounds on the number of arithmetic/Boolean operations that these algorithms require. We state several open questions including the existence of polynomial algorithms to resolve the above problems and availability of algorithms for addressing complete positivity over the rationals and over {0, 1} matrices.     

Palindromic quadratization and structure-preserving algorithm for palindromic matrix polynomials of even degree

In this paper, we propose a palindromic quadratization approach, transforming a palindromic matrix polynomial of even degree to a palindromic quadratic pencil. Based on the (S+ S^-1){(\mathcal{S}+ \mathcal{S}^{-1})} -transform and Patel’s algorithm, the structure-preserving algorithm can then be applied to solve the corresponding palindromic quadratic eigenvalue problem. Numerical experiments show that the relative residuals for eigenpairs of palindromic polynomial eigenvalue problems computed by palindromic quadratized eigenvalue problems are better than those via palindromic linearized eigenvalue problems or polyeig{{\texttt {polyeig}}} in MATLAB.     

Fast enclosure for the minimum norm least squares solution of the matrix equation AXB = C

Fast algorithms for enclosing the minimum norm least squares solution of the matrix equation AXB = C are proposed. To develop these algorithms, theories for obtaining error bounds of numerical solutions are established. The error bounds obtained by these algorithms are verified in the sense that all the possible rounding errors have been taken into account. Techniques for accelerating the enclosure and obtaining smaller error bounds are introduced. Numerical results show the properties of the proposed algorithms. Copyright © 2015 John Wiley & Sons, Ltd.     

Direct algorithms to solve the two‐sided obstacle problem for an M‐matrix

In this paper, two direct algorithms for solving the two‐sided obstacle problem with an M‐matrix are presented. The algorithms are well defined and have polynomial computational complexity. Copyright © 2006 John Wiley & Sons, Ltd.     

Orthogonal polynomial expansions for the matrix exponential

Many different algorithms have been suggested for computing the matrix exponential. In this paper, we put forward the idea of expanding in either Chebyshev, Legendre or Laguerre orthogonal polynomials. In order for these expansions to converge quickly, we cluster the eigenvalues into diagonal blocks and accelerate using shifting and scaling.     

Error-free transformations of matrix multiplication by using fast routines of matrix multiplication and its applications

This paper is concerned with accurate matrix multiplication in floating-point arithmetic. Recently, an accurate summation algorithm was developed by Rump et al. (SIAM J Sci Comput 31(1):189–224, 2008). The key technique of their method is a fast error-free splitting of floating-point numbers. Using this technique, we first develop an error-free transformation of a product of two floating-point matrices into a sum of floating-point matrices. Next, we partially apply this error-free transformation and develop an algorithm which aims to output an accurate approximation of the matrix product. In addition, an a priori error estimate is given. It is a characteristic of the proposed method that in terms of computation as well as in terms of memory consumption, the dominant part of our algorithm is constituted by ordinary floating-point matrix multiplications. The routine for matrix multiplication is highly optimized using BLAS, so that our algorithms show a good computational performance. Although our algorithms require a significant amount of working memory, they are significantly faster than ‘gemmx’ in XBLAS when all sizes of matrices are large enough to realize nearly peak performance of ‘gemm’. Numerical examples illustrate the efficiency of the proposed method.     

Numerical stability of biorthogonal wavelet transforms

Biorthogonal wavelets are essential tools for numerous practical applications. It is very important that wavelet transforms work numerically stable in floating point arithmetic. This paper presents new results on the worst-case analysis of roundoff errors occurring in floating point computation of periodic biorthogonal wavelet transforms, i.e. multilevel wavelet decompositions and reconstructions. Both of these wavelet algorithms can be realized by matrix–vector products with sparse structured matrices. It is shown that under certain conditions the wavelet algorithms can be remarkably stable. Numerous tests demonstrate the performance of the results.      

A heuristic for decomposing traffic matrices in TDMA satellite communication

A heuristic for decomposing traffic matrices in TDMA satellite communication. With the time-division multiple access (TDMA) technique in satellite communication the problem arises to decompose a givenn×n traffic matrix into a weighted sum of a small number of permutation matrices such that the sum of the weights becomes minimal. There are polynomial algorithms when the number of permutation matrices in a decomposition is allowed to be as large asn ². When the number of matrices is restricted ton, the problem is NP-hard. In this paper we propose a heuristic based on a scaling technique which for each number of allowed matrices in the range fromn ton ² allows to give a performance guarantee with respect to the total weight of the solution. As a subroutine we use new heuristic methods for decomposing a matrix of small integers into as few matrices as possible without exceeding the lower bound on the total weight. Computational results indicate that the method might also be practical.This work was supported by the Fonds zur Förderung der wissenschaftlichen Forschung, Project S32/01.     

Bezoutian矩阵的一致逼近形式

借助闭区间上的连续函数可以用Bernstein 多项式一致逼近这一事实,将多项式对所生成的经典Bezoutian 矩阵和Bernstein Bezoutian 矩阵推广到C [0,1]上函数对所对应的情形,给出了 Bezoutian 矩阵一致逼近形式的定义,并且得到如下结论：给出了经典 Bezoutian 矩阵的 Barnett 型分解公式和三角分解公式的一致逼近形式;提供了经典Bezoutian 矩阵和Bernstein Bezoutian 矩阵的一致逼近形式的两类算法;得到了上述两种矩阵的一致逼近形式中元素间的两个恒等关系式。最后,利用数值实例对恒等关系式进行验证,结果表明两类算法是有效的。     

Using FFT-based techniques in polynomial and matrix computations: recent advances and applicatons

Computations with univariate polynomials, like the evaluation of product, quotient, remainder, greatest common divisor, etc, are closely related to linear algebra computations performed with structured matrices having the Toeplitz-like or the Hankel-like structures.

The discrete Fourier transform, and the FFT algorithms for its computation, constitute a powerful tool for the design and analysis of fast algorithms for solving problems involving polynomials and structured matrices.

We recall the main correlations between polynomial and matrix computations and present two recent results in this field: in particular, we show how Fourier methods can speed up the solution of a wide class of problems arising in queueing theory where certain Markov chains, defined by infinite block Toeplitz matrices in generalized Hessenberg form, must be solved. Moreover, we present a new method for the numerical factorization of polynomials based on a matrix generalization of Koenig's theorem. This method, that relies on the evaluation/interpolation technique at the Fourier points, reduces the problem of polynomial factorization to the computation of the LU decomposition of a banded Toeplitz matrix with its rows and columns suitably permuted. Numerical experiments that show the effectiveness of our algorithms are presented     

On the regularity analysis of interpolatory Hermite subdivision schemes

It is well known that the critical Hölder regularity of a subdivision schemes can typically be expressed in terms of the joint-spectral radius (JSR) of two operators restricted to a common finite-dimensional invariant subspace. In this article, we investigate interpolatory Hermite subdivision schemes in dimension one and specifically those with optimal accuracy orders. The latter include as special cases the well-known Lagrange interpolatory subdivision schemes by Deslauriers and Dubuc. We first show how to express the critical Hölder regularity of such a scheme in terms of the joint-spectral radius of a matrix pair {F₀,F₁} given in a very explicit form. While the so-called finiteness conjecture for JSR is known to be not true in general, we conjecture that for such matrix pairs arising from Hermite interpolatory schemes of optimal accuracy orders a “strong finiteness conjecture” holds: ρ(F₀,F₁)=ρ(F₀)=ρ(F₁). We prove that this conjecture is a consequence of another conjectured property of Hermite interpolatory schemes which, in turn, is connected to a kind of positivity property of matrix polynomials. We also prove these conjectures in certain new cases using both time and frequency domain arguments; our study here strongly suggests the existence of a notion of “positive definiteness” for non-Hermitian matrices.     

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.     

Fast verification of solutions of matrix equations

Summary. In this paper, we are concerned with a matrix equation where A is an real matrix and x and b are n-vectors. Assume that an approximate solution is given together with an approximate LU decomposition. We will present fast algorithms for proving nonsingularity of A and for calculating rigorous error bounds for . The emphasis is on rigour of the bounds. The purpose of this paper is to propose different algorithms, the fastest with flops computational cost for the verification step, the same as for the LU decomposition. The presented algorithms exclusively use library routines for LU decomposition and for all other matrix and vector operations. Received June 16, 1999 / Revised version received January 25, 2001 / Published online June 20, 2001     

On the complexity of general matrix scaling and entropy minimization via the RAS algorithm

Given an n × m nonnegative matrix A = (a _ij ) and positive integral vectors and having a common one-norm h, the (r,c)-scaling problem is to obtain positive diagonal matrices X and Y, if they exist, such that XAY has row and column sums equal to r and c, respectively. The entropy minimization problem corresponding to A is to find an n × m matrix z = (z _ij ) having the same zero pattern as A, the sum of whose entries is a given number h, its row and column sums are within given integral vectors of lower and upper bounds, and such that the entropy function consisting of the sum of the terms z _ij ln (z _ij /a _ij ) is minimized. When the lower and upper bounds coincide, matrix scaling and entropy minimization are closely related. In this paper we present several complexity bounds for the -approximate (r,c)-scaling problem, polynomial in n,m,h, , and ln , where V and v are the largest and the smallest positive entries of A, respectively. These bounds, although not polynomial in , not only provide alternative complexities for the polynomial time algorithms, but could result in better overall complexities. In particular, our theoretical analysis supports the practicality of the well-known RAS algorithm. In our analysis we obtain bounds on the norm of scaling vectors which will be used in deriving not only some of the above complexities, but also a complexity for square nonnegative matrices having positive permanent. In particular, our results extend, nontrivially, many bounds for the doubly stochastic scaling of square nonnegative matrices previously given by Kalantari and Khachiyan to the case of general (r,c)-scaling. Finally, we study a more general entropy minimization problem where row and column sums are constrained to lie in prescribed intervals, and the sum of all entries is also prescribed. Balinski and Demange described an RAS type algorithm for its continuous version, but did not analyze its complexity. We show that this algorithm produces an -approximate solution within complexity polynomial in n, m, h, and .