Fast algorithms for hierarchically semiseparable matrices

Semiseparable matrices and many other rank‐structured matrices have been widely used in developing new fast matrix algorithms. In this paper, we generalize the hierarchically semiseparable (HSS) matrix representations and propose some fast algorithms for HSS matrices. We represent HSS matrices in terms of general binary HSS trees and use simplified postordering notation for HSS forms. Fast HSS algorithms including new HSS structure generation and HSS form Cholesky factorization are developed. Moreover, we provide a new linear complexity explicit ULV factorization algorithm for symmetric positive definite HSS matrices with a low‐rank property. The corresponding factors can be used to solve the HSS systems also in linear complexity. Numerical examples demonstrate the efficiency of the algorithms. All these algorithms have nice data locality. They are useful in developing fast‐structured numerical methods for large discretized PDEs (such as elliptic equations), integral equations, eigenvalue problems, etc. Some applications are shown. Copyright © 2009 John Wiley & Sons, Ltd.     

Blended kernel approximation in the ℋ︁‐matrix techniques

Several types of ??‐matrices were shown to provide a data‐sparse approximation of non‐local (integral) operators in FEM and BEM applications. The general construction is applied to the operators with asymptotically smooth kernel function provided that the Galerkin ansatz space has a hierarchical structure. The new class of ??‐matrices is based on the so‐called blended FE and polynomial approximations of the kernel function and leads to matrix blocks with a tensor‐product of block‐Toeplitz (block‐circulant) and rank‐k matrices. This requires the translation (rotation) invariance of the kernel combined with the corresponding tensor‐product grids. The approach allows for the fast evaluation of volume/boundary integral operators with possibly non‐smooth kernels defined on canonical domains/manifolds in the FEM/BEM applications. (Here and in the following, we call domains canonical if they are obtained by translation or rotation of one of their parts, e.g. parallelepiped, cylinder, sphere, etc.) In particular, we provide the error and complexity analysis for blended expansions to the Helmholtz kernel. Copyright © 2002 John Wiley & Sons, Ltd.     

A technique for constructing non‐embeddable quasi‐residual designs

We propose a technique for constructing two infinite families of non‐embeddable quasi‐residual designs as soon as one such design satisfying certain conditions exists. The main tools are generalized Hadamard matrices and balanced generalized weighing matrices. Starting with a specific non‐embeddable quasi‐residual 2‐(27,9,4) design, we construct for every positive integer m a non‐embeddable 2‐(3^m,3^m?1,(3^m?1?1)/2)‐design, and, if r_m=(3^m?1)/2 is a prime power, we construct for every positive integer n a non‐embeddable design. For each design in these families, a symmetric design with the corresponding parameters is known to exist. © 2002 Wiley Periodicals, Inc. J Combin Designs 10: 160–172, 2002; Published online in Wiley InterScience ( www.interscience.wiley.com ). DOI 10.1002/jcd.900     

Directional ‐matrix compression for high‐frequency problems

Standard numerical algorithms, such as the fast multipole method or ‐matrix schemes, rely on low‐rank approximations of the underlying kernel function. For high‐frequency problems, the ranks grow rapidly as the mesh is refined, and standard techniques are no longer attractive. Directional compression techniques solve this problem by using decompositions based on plane waves. Taking advantage of hierarchical relations between these waves' directions, an efficient approximation is obtained. This paper is dedicated to directional ‐matrices that employ local low‐rank approximations to handle directional representations efficiently. The key result is an algorithm that takes an arbitrary matrix and finds a quasi‐optimal approximation of this matrix as a directional ‐matrix using a prescribed block tree. The algorithm can reach any given accuracy, and the approximation requires only units of storage, where n is the matrix dimension, κ is the wave number, and k is the local rank. In particular, we have a complexity of if κ is constant and for high‐frequency problems characterized by κ²～n. Because the algorithm can be applied to arbitrary matrices, it can serve as the foundation of fast techniques for constructing preconditioners.     

Optimal (r,N)‐policy for discrete‐time Geo ∕ G ∕ 1 queue with different input rate and setup time

This paper presents a queue‐length analysis of Geo ∕ G ∕ 1 queue with ( r , N )‐policy and different input rate. Using a different method, the recursive expressions of queue‐length distribution at different epochs are obtained. Furthermore, some performance measures are also investigated. Finally, the Tabu search algorithm is used to search the joint optimum value of ( r , N ), which minimizes the state‐dependent operating cost. Copyright © 2014 John Wiley & Sons, Ltd.     

Compact representation of the full Broyden class of quasi‐Newton updates

In this paper, we present the compact representation for matrices belonging to the Broyden class of quasi‐Newton updates, where each update may be either rank one or rank two. This work extends previous results solely for the restricted Broyden class of rank‐two updates. In this article, it is not assumed that the same Broyden update is used in each iteration; rather, different members of the Broyden class may be used in each iteration. Numerical experiments suggest that a practical implementation of the compact representation is able to accurately represent matrices belonging to the Broyden class of updates. Furthermore, we demonstrate how to compute the compact representation for the inverse of these matrices and a practical algorithm for solving linear systems with members of the Broyden class of updates. We demonstrate through numerical experiments that the proposed linear solver is able to efficiently solve linear systems with members of the Broyden class of matrices with high accuracy. As an immediate consequence of this work, it is now possible to efficiently compute the eigenvalues of any limited‐memory member of the Broyden class of matrices, allowing for the computation of condition numbers and the ability to perform sensitivity analysis.     

A note on the complex semi‐definite matrix Procrustes problem

This note outlines an algorithm for solving the complex ‘matrix Procrustes problem’. This is a least‐squares approximation over the cone of positive semi‐definite Hermitian matrices, which has a number of applications in the areas of Optimization, Signal Processing and Control. The work generalizes the method of Allwright (SIAM J. Control Optim. 1988; 26 (3):537–556), who obtained a numerical solution to the real‐valued version of the problem. It is shown that, subject to an appropriate rank assumption, the complex problem can be formulated in a real setting using a matrix‐dilation technique, for which the method of Allwright is applicable. However, this transformation results in an over‐parametrization of the problem and, therefore, convergence to the optimal solution is slow. Here, an alternative algorithm is developed for solving the complex problem, which exploits fully the special structure of the dilated matrix. The advantages of the modified algorithm are demonstrated via a numerical example. Copyright © 2007 John Wiley & Sons, Ltd.     

A new efficient approach to the characterization of D‐stable matrices

The concept of D‐stability is relevant for stable square matrices of any order, especially when they appear in ordinary differential systems modeling physical problems. Indeed, D‐stability was treated from different points of view in the last 50 years, but the problem of characterization of a general D‐stable matrix was solved for low‐order matrices only (ie, up to order 4). Here, a new approach is proposed within the context of numerical linear algebra. Starting from a known necessary and sufficient condition, other simpler equivalent necessary and sufficient conditions for D‐stability are proved. Such conditions turn out to be computationally more appealing for symbolic software, as discussed in the reported examples. Therefore, a new symbolic method is proposed to characterize matrices of order greater than 4, and then it is used in some numerical examples, given in details.     

Quantized CP approximation and sparse tensor interpolation of function‐generated data

In this article, we consider the iterative schemes to compute the canonical polyadic (CP) approximation of quantized data generated by a function discretized on a large uniform grid in an interval on the real line. This paper continues the research on the quantics‐tensor train (QTT) method (“O(d log N)‐quantics approximation of N‐d tensors in high‐dimensional numerical modeling” in Constructive Approximation, 2011) developed for the tensor train (TT) approximation of the quantized images of function related data. In the QTT approach, the target vector of length 2^L is reshaped to a Lth‐order tensor with two entries in each mode (quantized representation) and then approximated by the QTT tensor including 2r²L parameters, where r is the maximal TT rank. In what follows, we consider the alternating least squares (ALS) iterative scheme to compute the rank‐r CP approximation of the quantized vectors, which requires only 2rL?2^L parameters for storage. In the earlier papers (“Tensors‐structured numerical methods in scientific computing: survey on recent advances” in Chemom Intell Lab Syst, 2012), such a representation was called Q_Can format, whereas in this paper, we abbreviate it as the QCP (quantized canonical polyadic) representation. We test the ALS algorithm to calculate the QCP approximation on various functions, and in all cases, we observed the exponential error decay in the QCP rank. The main idea for recovering a discretized function in the rank‐r QCP format using the reduced number of the functional samples, calculated only at O(2rL) grid points, is presented. The special version of the ALS scheme for solving the arising minimization problem is described. This approach can be viewed as the sparse QCP‐interpolation method that allows to recover all 2rL representation parameters of the rank‐r QCP tensor. Numerical examples show the efficiency of the QCP‐ALS‐type iteration and indicate the exponential convergence rate in r.     

Split‐type octonion matrix

The split and hyperbolic (countercomplex) octonions are eight‐dimensional nonassociative algebras over the real numbers, which are in the form , where e_m's have different properties for them. The main purpose of this paper is to define the split‐type octonion and its matrix whose inputs are split‐type octonions and give some properties for them by using the real quaternions, split, and hyperbolic (countercomplex) octonions. On the other hand, to make some definitions, we present some operations on the split‐type octonions. Also, we show that every split‐type octonions can be represented by 2 × 2 real quaternion matrix and 4 × 4 complex number matrix. The information about the determinants of these matrix representations is also given. Besides, the main features of split‐type octonion matrix concept are given by using properties of  real quaternion matrices. Then, 8n × 8nreal matrix representations of split‐type octonion matrices are shown, and some algebraic structures are examined. Additionally, we introduce real quaternion adjoint matrices of split‐type octonion matrices. Moreover, necessary and sufficient conditions and definitions are given for split‐type octonion matrices to be special split‐type octonion matrices. We describe some special split‐type octonion matrices. Finally, oct‐determinant of split‐type octonion matrices is defined. Definitive and understandable examples of all definitions, theorems, and conclusions were given for a better understanding of all these concepts.     

Automorphisms of higher‐dimensional Hadamard matrices

This article derives from first principles a definition of equivalence for higher‐dimensional Hadamard matrices and thereby a definition of the automorphism group for higher‐dimensional Hadamard matrices. Our procedure is quite general and could be applied to other kinds of designs for which there are no established definitions for equivalence or automorphism. Given a two‐dimensional Hadamard matrix H of order ν, there is a Product Construction which gives an order ν proper n‐dimensional Hadamard matrix P⁽ⁿ⁾(H). We apply our ideas to the matrices P⁽ⁿ⁾(H). We prove that there is a constant c > 1 such that any Hadamard matrix H of order ν > 2 gives rise via the Product Construction to cν inequivalent proper three‐dimensional Hadamard matrices of order ν. This corrects an erroneous assertion made in the literature that ”P⁽ⁿ⁾(H) is equivalent to “P⁽ⁿ⁾(H′) whenever H is equivalent to H′.” We also show how the automorphism group of P⁽ⁿ⁾(H) depends on the structure of the automorphism group of H. As an application of the above ideas, we determine the automorphism group of P⁽ⁿ⁾(Hk) when Hk is a Sylvester Hadamard matrix of order 2^k. For ν = 4, we exhibit three distinct families of inequivalent Product Construction matrices P⁽ⁿ⁾(H) where H is equivalent to H₂. These matrices each have large but non‐isomorphic automorphism groups. © 2008 Wiley Periodicals, Inc. J Combin Designs 16: 507–544, 2008     

On the solutions of the equation AXB = C under Toeplitz‐like and Hankel matrices constraint

In this paper, we are mainly concerned with 2 types of constrained matrix equation problems of the form AXB=C, the least squares problem and the optimal approximation problem, and we consider several constraint matrices, such as general Toeplitz matrices, upper triangular Toeplitz matrices, lower triangular Toeplitz matrices, symmetric Toeplitz matrices, and Hankel matrices. In the first problem, owing to the special structure of the constraint matrix , we construct special algorithms; necessary and sufficient conditions are obtained about the existence and uniqueness for the solutions. In the second problem, we use von Neumann alternating projection algorithm to obtain the solutions of problem. Then we give 2 numerical examples to demonstrate the effectiveness of the algorithms.     

When is a matrix unitary or Hermitian plus low rank?

Hermitian and unitary matrices are two representatives of the class of normal matrices whose full eigenvalue decomposition can be stably computed in quadratic computing complexity once the matrix has been reduced, for instance, to tridiagonal or Hessenberg form. Recently, fast and reliable eigensolvers dealing with low‐rank perturbations of unitary and Hermitian matrices have been proposed. These structured eigenvalue problems appear naturally when computing roots, via confederate linearizations, of polynomials expressed in, for example, the monomial or Chebyshev basis. Often, however, it is not known beforehand whether or not a matrix can be written as the sum of a Hermitian or unitary matrix plus a low‐rank perturbation. In this paper, we give necessary and sufficient conditions characterizing the class of Hermitian or unitary plus low‐rank matrices. The number of singular values deviating from 1 determines the rank of a perturbation to bring a matrix to unitary form. A similar condition holds for Hermitian matrices; the eigenvalues of the skew‐Hermitian part differing from 0 dictate the rank of the perturbation. We prove that these relations are linked via the Cayley transform. Then, based on these conditions, we identify the closest Hermitian or unitary plus rank k matrix to a given matrix A, in Frobenius and spectral norm, and give a formula for their distance from A. Finally, we present a practical iteration to detect the low‐rank perturbation. Numerical tests prove that this straightforward algorithm is effective.     

On some subclasses of P‐matrices

A matrix with positive row sums and all its off‐diagonal elements bounded above by their corresponding row averages is called a B‐matrix by J. M. Peña in References (SIAM J. Matrix Anal. Appl. 2001; 22 :1027–1037) and (Numer. Math. 2003; 95 :337–345). In this paper, it is generalized to more extended matrices—MB‐matrices, which is proved to be a subclass of the class of P‐matrices. Subsequently, we establish relationships between defined and some already known subclasses of P‐matrices (see, References SIAM J. Matrix Anal. Appl. 2000; 21 :67–78; Linear Algebra Appl. 2004; 393 :353–364; Linear Algebra Appl. 1995; 225 :117–123). As an application, some subclasses of P‐matrices are used to localize the real eigenvalues of a real matrix. Copyright © 2007 John Wiley & Sons, Ltd.     

Towards an optimal condition number of certain augmented Lagrangian‐type saddle‐point matrices

We present an analysis for minimizing the condition number of nonsingular parameter‐dependent 2 × 2 block‐structured saddle‐point matrices with a maximally rank‐deficient (1,1) block. The matrices arise from an augmented Lagrangian approach. Using quasidirect sums, we show that a decomposition akin to simultaneous diagonalization leads to an optimization based on the extremal nonzero eigenvalues and singular values of the associated block matrices. Bounds on the condition number of the parameter‐dependent matrix are obtained, and we demonstrate their tightness on some numerical examples. Copyright © 2016 John Wiley & Sons, Ltd.     

Two‐level discretization of the Navier‐Stokes equations with r‐Laplacian subgridscale viscosity

The r‐Laplacian has played an important role in the development of computationally efficient models for applications, such as numerical simulation of turbulent flows. In this article, we examine two‐level finite element approximation schemes applied to the Navier‐Stokes equations with r‐Laplacian subgridscale viscosity, where r is the order of the power‐law artificial viscosity term. In the two‐level algorithm, the solution to the fully nonlinear coarse mesh problem is utilized in a single‐step linear fine mesh problem. When modeling parameters are chosen appropriately, the error in the two‐level algorithm is comparable to the error in solving the fully nonlinear problem on the fine mesh. We provide rigorous numerical analysis of the two‐level approximation scheme and derive scalings which vary based on the coefficient r, coarse mesh size H, fine mesh size h, and filter radius δ. We also investigate the two‐level algorithm in several computational settings, including the 3D numerical simulation of flow past a backward‐facing step at Reynolds number Re = 5100. In all numerical tests, the two‐level algorithm was proven to achieve the same order of accuracy as the standard one‐level algorithm, at a fraction of the computational cost. © 2011 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2011     

The structured total least‐squares approach for non‐linearly structured matrices

In this paper, an extension of the structured total least‐squares (STLS) approach for non‐linearly structured matrices is presented in the so‐called ‘Riemannian singular value decomposition’ (RiSVD) framework. It is shown that this type of STLS problem can be solved by solving a set of Riemannian SVD equations. For small perturbations the problem can be reformulated into finding the smallest singular value and the corresponding right singular vector of this Riemannian SVD. A heuristic algorithm is proposed. Some examples of Vandermonde‐type matrices are used to demonstrate the improved accuracy of the obtained parameter estimator when compared to other methods such as least squares (LS) or total least squares (TLS). Copyright © 2002 John Wiley & Sons, Ltd.     

Computing a logarithm of a unitary matrix with general spectrum

We analyze an algorithm for computing a skew‐Hermitian logarithm of a unitary matrix and also skew‐Hermitian approximate logarithms for nearly unitary matrices. This algorithm is very easy to implement using standard software, and it works well even for unitary matrices with no spectral conditions assumed. Certain examples, with many eigenvalues near ? 1, lead to very non‐Hermitian output for other basic methods of calculating matrix logarithms. Altering the output of these algorithms to force skew‐Hermitian output creates accuracy issues, which are avoided by the considered algorithm. A modification is introduced to deal properly with the J‐skew‐symmetric unitary matrices. Applications to numerical studies of topological insulators in two symmetry classes are discussed. Copyright © 2014 John Wiley & Sons, Ltd.     

Accurate computations of matrices with bidiagonal decomposition using methods for totally positive matrices

A class of sign‐symmetric P‐matrices including all nonsingular totally positive matrices and their inverses as well as tridiagonal nonsingular H‐matrices is presented and analyzed. These matrices present a bidiagonal decomposition that can be used to obtain algorithms to compute with high relative accuracy their singular values, eigenvalues, inverses, or their LDU factorization. Copyright © 2012 John Wiley & Sons, Ltd.     

Parallel Newton two‐stage methods based on ILU factorizations for nonlinear systems

Parallel iterative algorithms based on the Newton method and on two of its variants, the Shamanskii method and the Chord method, for solving nonlinear systems are proposed. These algorithms are based on two‐stage multisplitting methods where incomplete LU factorizations are considered as a mean of constructing the inner splittings. Convergence properties of these parallel methods are studied for H‐matrices. Computational results of these methods on two parallel computing systems are discussed. The reported experiments show the effectiveness of these methods. Copyright © 2006 John Wiley & Sons, Ltd.