Square Matrix Research Articles

The problem of small and large matrices in parallel Matrix Multiplication

In this paper we discuss the case in which, using the generalized Cannon's algorithm, it is possible to reduce communications in matrix multiplication. We then apply reduction of communications to the case in which we have to multiply large matrices, in particular rectangular matrices. Two strategies are proposed to solve the problem of multiplying two large squared matrices. For the case in which we have to deal with small matrices, some methods are proposed to use the entire number of processors.

A square matrix is congruent to its transpose

For any matrix X let X ′ denote its transpose. We show that if A is an n by n matrix over a field K , then A and A ′ are congruent over K , i.e., P ′ AP = A ′ for some P ∈GL n ( K ).

The asymptotic variance of the univariate PLS estimator

In this short note, we derive an expression for the asymptotic covariance matrix of the univariate partial least squares (PLS) estimator. In contrast to M.C. Denham [J. Chemometrics 11 (1997) 39], who provided a locally linear approximation based on a recursive definition of the estimator, we derive a more compact expression for the asymptotic covariance matrix by combining a standard convergence result with matrix differential calculus, in particular the approach of J.R. Magnus and H. Neudecker [Matrix Differential Calculus with Applications in Statistics and Econometrics, revised ed., Wiley, Chichester, UK, 1991]. We also describe some theoretical and practical aspects of calculating the asymptotic covariance matrix, and illustrate its use on spectroscopic data.

Uniqueness of matrix square roots under a numerical range condition

For A∈M n( C), let W( A) denote the numerical range of A. It is shown that if W( A)∩(−∞,0)=∅, then A has a unique square root B∈M n( C) with W( B) in the closed right half plane.

Parallel multiplication of a vector by a kronecker product of matrices

Efficient parallel algorithms are derived for the multiplication of a vector by a Kronecker product of matrices. For N square matrices of order ni, i = 1, ..., N, it is well-known that the problem can be solved with a total of Πi=1N ni × Σi=1N ni floating point multiplications. The starting point of our investigations is a multidimensional stepwise scheme which respects the previous optimal bound. We give some sequential codes facing the memory space restriction and the sparsity context. Two parallel schedules depending on the choice between data communication and computation redundancies are proposed. In both cases, the number of processors is rather special. However, algorithms are quite efficient as illustrated by some experimental results.

Approximating the Logarithm of a Matrix to Specified Accuracy

The standard inverse scaling and squaring algorithm for computing the matrix logarithm begins by transforming the matrix to Schur triangular form in order to facilitate subsequent matrix square root and Pade approximation computations. A transformation-free form of this method that exploits incomplete Denman--Beavers square root iterations and aims for a specified accuracy (ignoring roundoff) is presented. The error introduced by using approximate square roots is accounted for by a novel splitting lemma for logarithms of matrix products. The number of square root stages and the degree of the final Pade approximation are chosen to minimize the computational work. This new method is attractive for high-performance computation since it uses only the basic building blocks of matrix multiplication, LU factorization and matrix inversion.

Uniqueness of matrix square roots and an application

Let A∈M n( C) . Let σ(A) denote the spectrum of A, and F(A) the field of values of A. It is shown that if σ(A)∩(−∞,0]=∅, then A has a unique square root B∈M n( C) with σ(B) in the open right (complex) half plane. This result and Lyapunov's theorem are then applied to prove that if F(A)∩(−∞,0]=∅, then A has a unique square root with positive definite Hermitian part. We will also answer affirmatively an open question about the existence of a real square root B∈M n( R) for A∈M n( R) with F(A)∩(−∞,0]=∅, where the field of values of B is in the open right half plane.

Successive matrix squaring algorithm for parallel computing the weighted generalized inverse AMN+

We derive a successive matrix squaring (SMS) algorithm to approximate the weighted generalized inverse, which can be expressed in the form of successive squaring of a composite matrix T. Given an m by n matrix A with m≈ n, we show that the weighted generalized inverse of A can be computed in parallel time ranging from O(log n) to O(log 2 n) provided that there are enough processors to support matrix multiplication in time O(log n).

Fluctuation analysis of stochastic gradient identification of polynomial Wiener systems

This correspondence presents analytical results and Monte Carlo simulations for the fluctuation behavior of a stochastic gradient adaptive identification scheme. This scheme identifies a polynomial Wiener system (linear FIR filter followed by a static polynomial nonlinearity) for noisy output observations. The analysis includes (1) bounds and a recursion for the misadjustment matrix and (2) algorithm mean square stability regions. A diagonal step-size matrix for the adaptive coefficients is introduced to speed up convergence. The theoretical predictions of the fluctuation analysis are supported by Monte Carlo simulations.

Visual encoding of patterns is subject to dual-task interference.

Two stimuli were presented at varying stimulus onset asynchronies (SOAs), with each stimulus associated with a distinct task. The first stimulus was a tone at one of either two or three frequencies. In two conditions, the task associated with a tone was either a speeded two-alternative discrimination (2AD), or a speeded three-alternative discrimination (3AD) based on the pitch of the tone. In a third condition, subjects were told to ignore the tone. The second stimulus was a briefly exposed study matrix of red and black squares followed by a mask. After a fixed delay, the mask was replaced by a test matrix that was either the same or different from the study matrix. The task associated with the matrices was to indicate, with no speed pressure, whether the study matrix and the test matrix were the same or different. Results from each speeded AD condition showed that subject's accuracy in the matrix task decreased as the SOA between the tone and the study matrix decreased. This effect was larger for the 3AD tone task than with a 2AD tone task. In addition, within each speeded AD condition, longer RTs in the tone task were associated with lower accuracy in the matrix task. None of these effects was evident when the subjects were told to ignore the tone. These results suggest that encoding visual information can be subject to significant capacity limitations imposed by cross-modal multitasking.

Successive matrix squaring algorithm for computing the Drazin inverse

This paper derives a successive matrix squaring (SMS) algorithm to approximate the Drazin inverse which can be expressed in the form of successive squaring of a composite matrix T . Given an n by n matrix A , the study shows that the Drazin inverse of A can be computed in parallel time ranging from O( log n) to O( log 2n) provided that there are enough processors to support matrix multiplication in time O( log n) . The SMS algorithm is generalized to higher-order schemes, where the composite matrix is repeatedly raised to an integer power l⩾2. This form of expression leads to a simplified notation compared to that of earlier methods, we argue that there is no obvious advantage in choosing l other than 2. Our derived error bound for the approximation of A D is new.

Extended Krylov Subspaces: Approximation of the Matrix Square Root and Related Functions

We introduce an economical Gram--Schmidt orthogonalization on the extended Krylov subspace originated by actions of a symmetric matrix and its inverse. An error bound for a family of problems arising from the elliptic method of lines is derived. The bound shows that, for the same approximation quality, the diagonal variant of the extended subspaces requires about the square root of the dimension of the standard Krylov subspaces using only positive or negative matrix powers. An example of an application to the solution of a 2.5-D elliptic problem attests to the computational efficiency of the method for large-scale problems.

A Bound for the Matrix Square Root with Application to Eigenvector Perturbation

Let $H$ be positive definite, let $\Delta H$ be Hermitian, and suppose that $\|H^{-1/2} (\Delta H) H^{-1/2} \| =1$. It is shown that the least constant $c_n$ such that \[ || (H + \eta \Delta H)^{1/2}H^{-1/2} - I || \leq c_n \eta + 0(\eta^2) \] for all $n \times n$ $H$ and $\Delta H$ grows like $\log n$. This fact has consequences in eigenvector perturbation theory.

Stable iterations for the matrix square root

Any matrix with no nonpositive real eigenvalues has a unique square root for which every eigenvalue lies in the open right half-plane. A link between the matrix sign function and this square root is exploited to derive both old and new iterations for the square root from iterations for the sign function. One new iteration is a quadratically convergent Schulz iteration based entirely on matrix multiplication; it converges only locally, but can be used to compute the square root of any nonsingular M-matrix. A new Pade iteration well suited to parallel implementation is also derived and its properties explained. Iterative methods for the matrix square root are notorious for suffering from numerical instability. It is shown that apparently innocuous algorithmic modifications to the Pade iteration can lead to instability, and a perturbation analysis is given to provide some explanation. Numerical experiments are included and advice is offered on the choice of iterative method for computing the matrix square root.

Weak localization and the transmission matrix.

We study the weak localization effect on the distribution of transmission matrix square eigenvalues in the framework of a fully microscopic approach. We identify the universal features of the distribution change and ascribe those to the lowest Cooperon mode. Details of the distribution are sensitive to an anomalously small magnetic field.

Efficient computation of the eigensystem realization algorithm

The computational efficiency of an algebraically equivalent alternative formulation of the eigensystem realization algorithm is investigated. This alternative can process large data sets much faster than the form of the eigensystem realization algorithm that is ordinarily used. This is primarily due to the partial eigenvalue decomposition of the inner or outer product of the Hankel matrix and its transpose, whichever is smaller, in place of the full singular value decomposition of the Hankel matrix itself or of the Hankel matrix square. By computing only the largest subset of the Hankel matrix product eigenvalues, low-order models can be realized from very large data sets. The primary benefit of this is improved convergence on the model pole locations for practical modal identification problems using hundreds of sensors to identify hundreds of modal resonances. This paper investigates the computational speed increase over the original algorithm on experimental data and assesses the magnitude of model errors introduced by squaring the Hankel matrix. In addition, key computational strategies are presented for efficient implementation of the algorithm on modern computer architectures.

Correction to Extending the trend vector: The trend matrix and sample-based partial least squares

Correction to Extending the trend vector: The trend matrix and sample-based partial least squares

Sample Partial Autocorrelation Function of a Multivariate Time Series

The choice of a matrix square root in order to define a correlation coefficient is crucial for the notion of partial autocorrelation function (PACF) for a multivariate time series. Here this topic is revisited and, introducing a new matrix link coefficient between two random vectors, a general framework for estimating the PACF is given. This leads to new autoregressive estimation methods based on sample estimators of the partial autocorrelation coefficients. Moreover, some generalizations of the scalar Burg′s technique fit in this framework making the comparison of all these methods easier.

Extending the trend vector: the trend matrix and sample-based partial least squares.

Trend vector analysis [Carhart, R.E. et al., J. Chem. Inf. Comput. Sci., 25 (1985) 64], in combination with topological descriptors such as atom pairs, has proved useful in drug discovery for ranking large collections of chemical compounds in order of predicted biological activity. The compounds with the highest predicted activities, upon being tested, often show a several-fold increase in the fraction of active compounds relative to a randomly selected set. A trend vector is simply the one-dimensional array of correlations between the biological activity of interest and a set of properties or 'descriptors' of compounds in a training set. This paper examines two methods for generalizing the trend vector to improve the predicted rank order. The trend matrix method finds the correlations between the residuals and the simultaneous occurrence of descriptors, which are stored in a two-dimensional analog of the trend vector. The SAMPLS method derives a linear model by partial least squares (PLS), using the 'sample-based' formulation of PLS [Bush, B.L. and Nachbar, R.B., J. Comput.-Aided Mol. Design, 7 (1993) 587] for efficiency in treating the large number of descriptors. PLS accumulates a predictive model as a sum of linear components. Expressed as a vector of prediction coefficients on properties, the first PLS component is proportional to the trend vector. Subsequent components adjust the model toward full least squares. For both methods the residuals decrease, while the risk of overfitting the training set increases. We therefore also describe statistical checks to prevent overfitting. These methods are applied to two data sets, a small homologous series of disubstituted piperidines, tested on the dopamine receptor, and a large set of diverse chemical structures, some of which are active at the muscarinic receptor. Each data set is split into a training set and a test set, and the activities in the test set are predicted from a fit on the training set. Both the trend matrix and the SAMPLS approach improve the predictions over the simple trend vector. The SAMPLS approach is superior to the trend matrix in that it requires much less storage and CPU time. It also provides a useful set of axes for visualizing properties of the compounds. We describe a randomization method to determine the optimum number of PLS components that is very much faster for large training sets than leave-one-out cross-validation.

Generalised matrix inversion and rank computation by successive matrix powering

Computation of the generalised inverse A + and rank of an arbitrary (including singular and rectangular) matrix A has many applications. This paper derives an iterative scheme to approximate the generalised inverse which can be expressed in the form of successive squaring of a composite matrix T . Given an m by n matrix A with m≈n , we show that the generalised inverse of A can be computed in parallel time ranging from O (log n ) to O (log 2 n ), similar to previous methods. The rank of matrix A is obtained along with the generalised inverse. The successive matrix squaring algorithm is generalised to higher-order schemes, where the composite matrix is repeatedly raised to an integer power l >2. This form of expression leads to a simplified notation compared with that of earlier methods, and helps to clarify the relationship between l , the order of the iterative scheme and K , the number of iterations. In particular, the accuracy achieved in approximating A + is a function only of the magnitude of l K and does not depend on the particular values chosen for l and K ; we argue that there is no obvious advantage in choosing l other than 2. Our derived error bound for the approximation to A + is tighter than that previously established. The same bound applies to the rank. Numerical experiments with different test matrices (square, rectangular, complex, singular, etc.) illustrate the method. They further demonstrate that our tighter error bound provides a useful guide to the number of iterations required. In the examples given, the specified accuracy was achieved after the calculated number of iterations, but no earlier. Implementation studies on a general-purpose parallel machine (CM-5) demonstrated a smaller than expected penalty for direct squaring of matrix T in comparison with multiplication of its component block matrices. For special-purpose VLSI architectures, the simple structure of the matrix squaring algorithm leads to a straightforward parallel implementation with low communication overheads.