Matrix Methods: Theory, Algorithms and Applications

The following sections are included:

• PREFACE

• Contents

The following sections are included:

• Introduction

• Preliminaries

• Matrices over the binary Boolean semiring

• Matrices over antinegative semirings without zero divisors

• References

Since quaternions have isomorphic representations in matrix form we investigate various well known matrix decompositions for quaternions.

Stability of a linear autonomous non-conservative system in the presence of potential, gyroscopic, dissipative, and non-conservative positional forces is studied. The cases when the non-conservative system is close either to a gyroscopic system or to a circulatory one, are examined. It is known that marginal stability of gyroscopic and circulatory systems can be destroyed or improved up to asymptotic stability due to action of small non-conservative positional and velocity-dependent forces. We show that in both cases the boundary of the asymptotic stability domain of the perturbed system possesses singularities such as “Dihedral angle”, “Break of an edge” and “Whitney's umbrella” that govern stabilization and destabilization as well as are responsible for the imperfect merging of modes. Sensitivity analysis of the critical parameters is performed with the use of the perturbation theory for eigenvalues and eigenvectors of non-self-adjoint operators. In case of two degrees of freedom, stability boundary is found in terms of the invariants of matrices of the system. Bifurcation of the stability domain due to change of the structure of the damping matrix is described. As a mechanical example, the Hauger gyropendulum is analyzed in detail; an instability mechanism in a general mechanical system with two degrees of freedom, which originates after discretization of models of a rotating disc in frictional contact and possesses the spectral mesh in the plane ‘frequency’ versus ‘angular velocity’, is analytically described and its role in the excitation of vibrations in the squealing disc brake and in the singing wine glass is discussed.

For each square complex matrix, V. I. Arnold constructed a normal form with the minimal number of parameters to which a family of all matrices B that are close enough to this matrix can be reduced by similarity transformations that smoothly depend on the entries of B. Analogous normal forms were also constructed for families of complex matrix pencils by A. Edelman, E. Elmroth, and B. Kågström, and contragredient matrix pencils (i.e., of matrix pairs up to transformations (A, B) ↦ (S^-1AR, R^-1BS)) by M. I. Garcia-Planas and V. V. Sergeichuk. In this paper we give other normal forms for families of matrices, matrix pencils, and contragredient matrix pencils; our normal forms are block triangular.

In this paper we present some results of Schein rank of Boolean matrices. A notion of the intersection number of a bipartite graph is defined and its applications to Schein rank of Boolean matrices are derived. We discuss minimal and maximal matrices of given Schein rank, the number of m × n Boolean matrices with given Schein rank. The Schein ranks of some m × n Boolean matrices are determined. In the last section, we give some further result concerning the Schein rank of Boolean matrices.

We consider matrices over a Brouwerian lattice. The linear span of columns of a matrix A form a semilattice. We call it a column semilattice for A. The questions are: when column semilattice is a lattice, when column semilattice is a distributive lattice, and what formulas can be obtained for the meet and the join operations? We prove that for any lattice matrix A, the column semilattice is a lattice. We also obtain formulas for the meet and the join operations. If A is an idempotent or A is a regular matrix, then the column semilattice is a distributive lattice.

We also consider invariant eigenvectors of a square matrix A over a Brouwerian lattice. It is proved that all A-invariant eigenvectors form a distributive lattice and the simple formulas for the meet and the join operations are obtained.

Let be a field and let A be a finite-dimensional -algebra. We define the length of a finite generating set of this algebra as the smallest number k such that words of the length not greater than k generate A as a vector space, and the length of the algebra is the maximum of lengths of its generating sets. In this paper we study the connection between the length of an algebra and the lengths of its subalgebras. It turns out that the length of an algebra can be smaller than the length of its subalgebra. To investigate, how different the length of an algebra and the length of its subalgebra can be, we evaluate the difference and the ratio of the lengths of an algebra and its subalgebra for several representative families of algebras. Also we give examples of length computation of two and three block upper triangular matrix algebras.

The objective of this paper is to consider a class of singular nonsymmetric matrices with integer spectrum. The class comprises generalized triangular matrices with diagonal elements obtained by summing the elements of the corresponding column. If the size of a matrix belonging to the class equals n × n, the spectrum of the matrix is given by the sequence of distinct nonnegative integers up to n − 1, irrespective of the elements of the matrix. Right and left eigenvectors are obtained. Moreover, several interesting relations are presented, including factorizations via triangular matrices.

We show that for n × n nonsingular matrices A₁, A₂, …, A_k (k ≥ 2) over a commutative principal ideal domain R with relatively coprime determinants there exist invertible n × n matrices U, V₁, …, V_k over R such that UA_iV_i = S_Ai are the Smith normal forms of the matrices A_i.

We survey theoretical properties and algorithms concerning the problem of solving a nonsymmetric algebraic Riccati equation, and we report on some known methods and new algorithmic advances. In particular, some results on the number of positive solutions are proved and a careful convergence analysis of Newton's iteration is carried out in the cases of interest where some singularity conditions are encountered. From this analysis we determine initial approximations which still guarantee the quadratic convergence.

A new version of the generalized conjugate direction (GCD) method for nonsymmetric linear algebraic systems is proposed which is oriented on large and ill-conditioned sets of equations. In distinction from the known Krylov subspace methods for unsymmetrical matrices, the method uses explicitly computed A-conjugate (in generalized sense) vectors, along with an orthogonal set of residuals obtained in the Arnoldi orthogonalization process. Employing entire sequences of orthonormal basis vectors in the Krylov subspaces, similarly to GMRES and FOM, ensures high stability of the method. But instead of solution of a linear set of equations with a Hessenberg matrix in each iteration for determining the step we use A-conjugate vectors and some simple recurrence formulas. The performance of the proposed algorithm is illustrated by the results of extensive numerical experiments with large-scale ill-conditioned linear systems and by comparison with the known efficient algorithms.

We answer a question motivated by our study of the normal Hankel problem, i.e., the problem of describing normal Hankel matrices. It was shown previously that new solutions to this problem can only be found among the so-called (φ, ψ)-circulants. The latter can be described by a system of equations with respect to the real and imaginary parts of their entries. Since the equations are quadratic, it is not at all clear whether this system admits real solutions unless n (the order of the matrix) is three or four (these cases were solved in an earlier publication of the authors). In this note, we construct a class of normal Hankel matrices of an arbitrary order n ≥ 5 that are (φ, ψ)-circulants for appropriately chosen values of φ and ψ.

The abrupt boundary truncation of an image introduces artifacts in the restored image. For large image restoration with shift-invariant blurring, it is advisable to use Fast Fourier Transform (FFT)-based procedures for reducing the computational effort. In this direction several techniques manipulate the observed image at the boundary or make some assumptions on the boundary of the true image, in such a way that FFT-based algorithms can be used. We compare the use of reflection with that of anti-reflection, in connection with the choice of the boundary conditions or for extending the observed image, both theoretically and numerically. Furthermore, we combine the two proposals. More precisely we apply anti-reflection, followed by reflection if necessary, to the observed image and we observe that the resulting restoration quality is increased with respect to the case of plain reflection.

Jim Wilkinson discovered that the computation of zeros of polynomials is ill conditioned when the polynomial is given by its coefficients. For many problems we need to compute zeros of polynomials, but we do not necessarily need to represent the polynomial with its coefficients. We develop algorithms that avoid the coefficients. They turn out to be stable, however, the drawback is often heavily increased computational effort. Modern processors on the other hand are mostly idle and wait for crunching numbers so it may pay to accept more computations in order to increase stability and also to exploit parallelism. We apply the method for nonlinear eigenValue problems.

Pseudoskeleton approximation and some other problems require the knowledge of sufficiently well-conditioned submatrix in a largescale matrix. The quality of a submatrix can be measured by modulus of its determinant, also known as volume. In this paper we discuss a search algorithm for the maximum-volume submatrix which already proved to be useful in several matrix and tensor approximation algorithms. We investigate the behavior of this algorithm on random matrices and present some of its applications, including maximization of a bivariate functional.

The acceleration of the original projective iterative methods of multiplicative or additive type for solving systems of linear algebraic equations (SLAEs) by means of conjugate direction approaches is considered. The orthogonal and varitional properties of the preconditioned conjugate gradient, conjugate residual and semi-conjugate residual algorithms, as well as estimations of the number of iterations, are presented. Similar results were obtained for the dynamically preconditioned iterative process in Krylov subspaces. Application of discussed techniques for domain decomposition, Kaczmarz, and Cimmino methods is proposed.

For any given n-by-n matrix A_n, a specific circulant preconditioner t_F(A_n) proposed by E. Tyrtyshnikov [SIAM J. Matrix Anal. Appl., Vol. 13 (1992), pp. 459–473] is defined to be the solution of

over all n-by-n nonsingular circulant matrices C_n. The t_F(A_n), called the superoptimal circulant preconditioner, has been proved to be a good preconditioner for a large class of structured systems including some ill-conditioned problems from image processing. In this paper, we study this preconditioner from an operator viewpoint. We will give some relationships between the optimal preconditioner (operator) proposed by T. Chan [SIAM J. Sci. Statist. Comput., Vol. 9 (1988), pp. 766–771] and superoptimal preconditioner (operator).

A theoretical justification is found for several standard techniques related to ILU preconditioning, such as pre-scaling and pivot modification, with implications for practical implementation. An improved estimate for the reduction of the GMRES residual is obtained within the general framework of two-stage preconditioning. In particular, an estimate in terms of a conditioning measure of the scaled coefficient matrix and the Frobenius norm of the scaled ILU residual is presented.

In this paper, we re-investigate the resolution of Toeplitz systems Tu = g, from a new point of view, by correlating the solution of such problems with syzygies of polynomials or moving lines. We show an explicit connection between the generators of a Toeplitz matrix and the generators of the corresponding module of syzygies. We show that this module is generated by two elements of degree n and the solution of Tu = g can be reinterpreted as the remainder of an explicit vector depending on g, by these two generators.

This approach extends naturally to multivariate problems and we describe for Toeplitz-block-Toeplitz matrices, the structure of the corresponding generators.

We present concepts of data-sparse tensor approximations to the functions and operators arising in many-particle models of quantum chemistry. Our approach is based on the systematic use of structured tensor-product representations where the low-dimensional components are represented in hierarchical or wavelet based matrix formats. The modern methods of tensor-product approximation in higher dimensions are discussed with the focus on analytically based approaches. We give numerical illustrations which confirm the efficiency of tensor decomposition techniques in electronic structure calculations.

The transformation of the discrete nonlinear oscillatory system separating the Fourier coefficients and the harmonics of a linear oscillator system is obtained. The numerical experiments show that an initial impulse of the nonlinear system is localized in those harmonics whose numbers are divisible by the numbers of nonzero initial harmonics. Obtained separation of the variables is connected with symmetric Toeplitz and persymmetric Hankel matrices.

We accelerate multipoint polynomial evaluation by reducing the problem to structured matrix computation and transforming the resulting matrix structure.

We begin with specifying a class of matrices for which Gaussian elimination with partial pivoting fails and then observe that both rook and complete pivoting easily handle these matrices. We display the results of testing partial, rook and complete pivoting for this and other classes of matrices. Our tests confirm that rook pivoting is an inexpensive but solid backup wherever partial pivoting fails.

We first cover Newton's iteration for generalized matrix inversion, its ameliorations, recursive compression of its iterates in the case of structured inputs, some techniques of continuation via factorization, and extension to splitting the Singular Value Decomposition. We combine the latter extension with our recent fast algorithms for the null space bases (prompted by our progress in randomized preconditioning). We applied these combinations to compute the respective spaces of singular vectors and to arrive at divide-and-conquer algorithms for matrix inversion and computing determinants. Our techniques promise to be effective for computing other matrix functions in the case of ill conditioned inputs.

The paper analyzes and compares some spectral filtering methods as truncated singular/eigen-value decompositions and Tikhonov/Reblurring regularizations in the case of the recently proposed Reflective [18] and Anti-Reflective [21] boundary conditions. We give numerical evidence to the fact that spectral decompositions (SDs) provide a good image restoration quality and this is true in particular for the Anti-Reflective SD, despite the loss of orthogonality in the associated transform. The related computational cost is comparable with previously known spectral decompositions, and results substantially lower than the singular value decomposition. The model extension to the cross-channel blurring phenomenon of color images is also considered and the related spectral filtering methods are suitably adapted.

Polynomial matrices with positive semidefinite coefficients C_i are studied. If C₀ is positive definite and ∑ C_i = I then all characteristic values of G(z) are in the closed unit disc and those lying on the unit circle are m-th roots of unity having linear elementary divisors. The result yields a stability and convergence criterion for a system of difference equations.

We consider a boundary value problem whose generalized statement is formulated as a mixed variational inequality in a Hilbert space. The operator of this variational inequality is a sum of several inversely strongly monotone operators (which are not necessarily potential operators). The functional occurring in this variational inequality is also a sum of several lower semi-continuous convex proper functionals. For solving of the considered variational inequality a decomposition iterative method is offered. The suggested method does not require the inversion of original operators. The convergence of this method is investigated.

A class of multilevel algorithms for partitioning of a sparse matrix prior to parallel solution of a system of linear equations is described. This matrix partitioning problem can be described in terms of a graph partitioning problem which is known to be NP-hard, so several heuristics for its solution have been proposed in the past decades. For this purpose we use the multilevel algorithm proposed by B. Hendrickson and R. Leland [2] and further developed by G. Karypis and V. Kumar [3]. This algorithm is very efficient and tends to produce high quality partitioning for a wide range of matrices arising in many practical applications.

Singular Spectrum Analysis is a nonparametric method, which allows one to solve problems like decomposition of a time series into a sum of interpretable components, extraction of periodic components, noise removal and others. In this paper, the algorithm and theory of the SSA method are extended to analyse two-dimensional arrays (e.g. images). The 2D-SSA algorithm based on the SVD of a Hankel-block-Hankel matrix is introduced. Another formulation of the algorithm by means of Kronecker-product SVD is presented. Basic SSA notions such as separability are considered. Results on ranks of Hankel-block-Hankel matrices generated by exponential, sine-wave and polynomial 2D-arrays are obtained. An example of 2D-SSA application is presented.

A new approach for solution of the boundary value problems for wide class of elliptic partial differential equations of mathematical physics is proposed. This class includes the Laplace, Poisson, and Helmholtz equations. The approach discovered by author is based on the Local Ray Principle and leads to new General Ray (GR) method, which presents the solution of the Dirichlet boundary problems by explicit analytical formulas that include the direct and inverse Radon transform. GR-method is realized by fast algorithms and MATLAB software, whose quality is demonstrated by numerical experiments.

A new efficient multigrid algorithm is proposed for solving parabolic equations. It is similar to implicit schemes by stability and accuracy, but the computational complexity is substantially reduced at each time step. Stability and accuracy of the proposed two-grid algorithm are analyzed theoretically for one- and two-dimensional heat diffusion equations. Good accuracy is demonstrated on model problems for one- and two-dimensional heat diffusion equations, including those with thermal conductivity defined as a discontinuous function of coordinates.

A new finite volume scheme for 3D diffusion problems with heterogeneous full diffusion tensor is considered. The discretization uses nonlinear two-point flux approximation on unstructured tetrahedral grids. Monotonicity of the linearized operator allows us to guarantee non-negativity of the discrete solution.

We consider two-dimensional integro-differential equation for currents in thin superconducting films. The integral operator of this equation is hypersingular operator with kernel decaying as 1/R³. For numerical solution Galerkin Finite Element Method (FEM) on triangular mesh with linear elements is used. It results in dense FEM matrix of large dimension. As the kernel is quickly decaying then off-diagonal elements of FEM matrix are small. We investigate simple sparsification approach based on dropping small entries of FEM matrix. The conclusion is that it allows to reduce to some extent memory requirements. Nevertheless for problems with large number of mesh points more complicated techniques as one of hierarchical matrices algorithms should be considered.

The problem of the stationary magnetic field modelling in the presence of an ideal-conductive surface is reduced to the integro-differential equation of the first kind. The analysis of the equation is carried out by the variational method. The novel software package has been created for its solving. The examples of its usage are represented.

For RCLM networks we present a novel algebraic spectral model order reduction algorithm equipped with efficient tools for preserving the passivity. For RC networks our approach is similar to the well-known spectral reduction technique PACT (pole analysis via congruence transformations). The accuracy and reduction ratio of resulting reduced-order models are demonstrated with several industrial examples.

New smoothers resulted from a special class of triangular skew-symmetric splitting iteration methods for the multigrid methods were used to solve the systems of linear equations with strongly non-symmetric coefficient matrices, which may be produced by the central-difference approximation of a convection-diffusion equation with dominant convection.

The problem of eddy currents computation on the conducting surface is considered. It is reduced to the integral equations. The existence, uniqueness and stability of their solutions is proved. The numerical method for the mentioned problem is described. The example of computing is given.

The general vector boundary-value problem for the integro-differential kinetic Boltzmann equation, which describes the polarized radiation transfer in a heterogeneous plane layer with a horizontally non-homogeneous and anisotropically reflective boundary, cannot be solved by the finite-difference methods. A mathematical model is proposed and justified that gives an asymptotically exact solution of the boundary-value problem in the class of the functions of slow growth. The new model is constructed by the influence function method and is efficient for algorithms using parallel computing.

The method of the solving of ill-posed problems turned into the solving of arbitrary systems of linear algebraic equations is considered. This method is based on the reduction of an arbitrary (in general, inconsistent) linear system to an equivalent consistent augmented system with a symmetric matrix. The problem of choosing of regularization parameter is considered.

The following sections are included:

• Author Index