Fast Array Algorithms for Structured Matrices

Computer Solution of Large Linear Systems

Triangular factorization of general N × N Hermitian matrices requires O(N 3) operations. We present a procedure that factors certain structured matrices in O(N 2) operations. The family of structured matrices considered in this paper includes Toeplitz and Hankel matrices, as well as a variety of ‘Toeplitz-like’ and ‘Hankel-like’ matrices (e.g. inverses, sums and products of Toeplitz and Hankel matrices).

In this article, the fast algorithms for the discrete Tchebichef transform (DTT) are proposed for input sequences of lengths in the range from 3 to 8. At present, DTT is widely applied in signal processing, image compression, and video coding. The review of the articles related to fast DTT algorithms has shown that such algorithms are mainly developed for input signal lengths 4 and 8. However, several problems exist for which signal and image processing with different apertures is required. To avoid this shortcoming, the structural approach and a sparse matrix factorization are applied in this paper to develop fast real DTT algorithms for short-length input signals. According to the structural approach, the rows and columns of the transform matrix are rearranged, possibly by changing the signs of some rows or columns. Next, the matched submatrix templates are extracted from the matrix structure and decomposed into a matrix product to construct the factorization of an initial matrix. A sparse matrix factorization assumes that the butterfly architecture can be extracted from the transform matrix. Combining the structural approach with a sparse matrix factorization, we obtained the matrix representation with reduced computational complexity. Based on the obtained matrix representation, the fast algorithms were developed for the real DTT via the data flow graphs. The fast algorithms for integer DTT can be easily obtained using the constructed data flow graphs. To confirm the correctness of the designed algorithms, the MATLAB R2023b software was applied. The constructed factorizations of the real DTT matrices reduce the number of multiplication operations by 78% on average compared to the direct matrix-vector product at signal lengths in the range from 3 to 8. The number of additions decreased by 5% on average within the same signal length range.

Rank structures provide an opportunity to develop new efficient numerical methods for practical problems, when the off-diagonal blocks of certain dense intermediate matrices have small (numerical) ranks. In this work, we present a framework of structured direct factorizations for general sparse matrices, including discretized PDEs on general meshes, based on the multifrontal method and hierarchically semiseparable (HSS) matrices. We prove the idea of replacing certain complex structured operations by fast simple ones performed on compact reduced matrix forms. Such forms result from the hierarchical factorization of a tree-structured HSS matrix in a ULV-type scheme, so that the tree structure is reduced into a single node, the root of the original tree. This idea is shown to be very useful in the partial ULV factorization of an HSS matrix (for quickly computing Schur complements) as well as the solution stage. These techniques are then built into the multifrontal method for sparse factorizations after nested dissection, so as to convert the intermediate dense factorizations into fast structured ones. This method keeps certain Schur complements dense so as to avoid complicated data assembly, and is much simpler and more general than some existing methods. In particular, if the matrix arises from the discretization of certain PDEs, the factorization costs roughly $O(n)$ flops in two dimensions, and roughly $O(n^{4/3})$ flops or less in three dimensions. The solution cost and memory are nearly $O(n)$ in both cases. These counts are obtained with an idea of rank relaxation, so that this method is more generally applicable to problems where the intermediate off-diagonal ranks are not small. We demonstrate the performance of the method with two- and three-dimensional discretized equations, as well as various examples from a sparse matrix collection. The ideas here are also useful in future developments of fast structured solvers.

To accurately estimate multi-section traffic states of freeway, a new multi-section error state equations are built up. In the new model, error propagations of both the traffic density and the average velocity are considered and the extended Kalman filter is used to estimate all traffic states. To avoid a series of high-dimension problems, a modified weighted Gram-Schmidt orthogonal U-D factorization method is used for the time update and measurement update of the extended Kalman filter to get high numerical stability and computational efficiency. Considered the structure of system matrix, block matrix is used in U-D factorization algorithm. Results of simulation to 100 section states of freeway and actual application show that the new method can be efficiently used to estimation and prediction freeway traffic flows.

The fast algorithms of discrete sinusoidal transform of the fourth type (DST-IV) for small-length input data in the range of lengths from 2 to 9 are developed. Fast algorithms for short input data sequences are subsequently used as building blocks for designing fast algorithms of large-sized discrete transforms. Applying the fast DST-IV algorithms for small-size block processing can reduce overall system complexity and delay, allowing detailed signal processing. As a result of the literature review, two main approaches to developing fast discrete sine transform (DST) algorithms were identified, namely, the polynomial algebraic approach and the matrix factorization approach. In the paper, the last approach is exploited. A matrix-vector product expression of the DST-IV is the starting point for designing the fast algorithms. Then based on the repetition and arranging of the matrix elements, the factorisation of the matrices of coefficients of DST-IV is produced to reduce computational complexity. The correctness of the obtained algorithmic solutions was justified theoretically using a strict mathematical background of each of them. The elaborated algorithms were then further tested using MATLAB R2023b software to finally confirm their performance. The resulting factorizations of the DST-IV matrices reduce the number of multiplications by 63% but increase the number of additions by 8% on average in the range of signal sample numbers from 3 to 9. It has been observed that for even-length input sequences, the reduction in the number of multiplications is not as significant as for odd-length sequences. For some other well-known discrete trigonometric transforms (discrete Fourier transform, discrete Hartley transform) the opposite situation holds. The proposed DST-IV fast algorithms do not limit the length of the input data sequence to powers of two or three. The data flow graphs constructed for the proposed algorithms reveal their modular space-time structure suitable for VLSI implementation. Keywords: discrete sine transform; matrix factorization; fast algorithms; computational complexity; digital signal processing.

Use of the stochastic Galerkin finite element methods leads to large systems\nof linear equations obtained by the discretization of tensor product solution\nspaces along their spatial and stochastic dimensions. These systems are\ntypically solved iteratively by a Krylov subspace method. We propose a\npreconditioner which takes an advantage of the recursive hierarchy in the\nstructure of the global matrices. In particular, the matrices posses a\nrecursive hierarchical two-by-two structure, with one of the submatrices block\ndiagonal. Each one of the diagonal blocks in this submatrix is closely related\nto the deterministic mean-value problem, and the action of its inverse is in\nthe implementation approximated by inner loops of Krylov iterations. Thus our\nhierarchical Schur complement preconditioner combines, on each level in the\napproximation of the hierarchical structure of the global matrix, the idea of\nSchur complement with loops for a number of mutually independent inner Krylov\niterations, and several matrix-vector multiplications for the off-diagonal\nblocks. Neither the global matrix, nor the matrix of the preconditioner need to\nbe formed explicitly. The ingredients include only the number of stiffness\nmatrices from the truncated Karhunen-Lo\\`{e}ve expansion and a good\npreconditioned for the mean-value deterministic problem. We provide a condition\nnumber bound for a model elliptic problem and the performance of the method is\nillustrated by numerical experiments.\n

Efficient parallel factorization and solution of structured and unstructured linear systems

Some generalized inverses of partition matrix and quotient identity of generalized Schur complement

Age-at-onset is one of the critical traits in cohort studies of age-related diseases. Large-scale genome-wide association studies (GWAS) of age-at-onset traits can provide more insights into genetic effects on disease progression and transitions between stages. Moreover, proportional hazards (or Cox) regression models can achieve higher statistical power in a cohort study than a case-control trait using logistic regression. Although mixed-effects models are widely used in GWAS to correct for sample dependence, application of Cox mixed-effects models (CMEMs) to large-scale GWAS is so far hindered by intractable computational cost. In this work, we propose COXMEG, an efficient R package for conducting GWAS of age-at-onset traits using CMEMs. COXMEG introduces fast estimation algorithms for general sparse relatedness matrices including, but not limited to, block-diagonal pedigree-based matrices. COXMEG also introduces a fast and powerful score test for dense relatedness matrices, accounting for both population stratification and family structure. In addition, COXMEG generalizes existing algorithms to support positive semidefinite relatedness matrices, which are common in twin and family studies. Our simulation studies suggest that COXMEG, depending on the structure of the relatedness matrix, is orders of magnitude computationally more efficient than coxme and coxph with frailty for GWAS. We found that using sparse approximation of relatedness matrices yielded highly comparable results in controlling false-positive rate and retaining statistical power for an ethnically homogeneous family-based sample. By applying COXMEG to a study of Alzheimer's disease (AD) with a Late-Onset Alzheimer's Disease Family Study from the National Institute on Aging sample comprising 3456 non-Hispanic whites and 287 African Americans, we identified the APOE ε4 variant with strong statistical power (P = 1e-101), far more significant than that reported in a previous study using a transformed variable and a marginal Cox model. Furthermore, we identified novel SNP rs36051450 (P = 2e-9) near GRAMD1B, the minor allele of which significantly reduced the hazards of AD in both genders. These results demonstrated that COXMEG greatly facilitates the application of CMEMs in GWAS of age-at-onset traits.

Let A and B be n-square positive definite matrices. Denote the Hadamard product of A and B by A∘ B. The main results of the paper are: 1.For any matrices C and D of size m × n and 2.Let A/α be the Schur complement of A(α) in A. Then . Some other matrix inequalities of Schur complements and Hadamard products of positive definite matrices are also presented.

This paper provides a detailed analysis that shows how to stabilize the {\em generalized} Schur algorithm, which is a fast procedure for the Cholesky factorization of positive-definite structured matrices $R$ that satisfy displacement equations of the form $R-FRF^T=GJG^T$, where $J$ is a $2\times 2$ signature matrix, $F$ is a stable lower-triangular matrix, and $G$ is a generator matrix. In particular, two new schemes for carrying out the required hyperbolic rotations are introduced and special care is taken to ensure that the entries of a Blaschke matrix are computed to high relative accuracy. Also, a condition on the smallest eigenvalue of the matrix, along with several computational enhancements, is introduced in order to avoid possible breakdowns of the algorithm by assuring the positive-definiteness of the successive Schur complements. We use a perturbation analysis to indicate the best accuracy that can be expected from {\em any} finite-precision algorithm that uses the generator matrix as the input data. We then show that the modified Schur algorithm proposed in this work essentially achieves this bound when coupled with a scheme to control the generator growth. The analysis further clarifies when pivoting strategies may be helpful and includes illustrative numerical examples. For all practical purposes, the major conclusion of the analysis is that the modified Schur algorithm is backward stable for a large class of structured matrices.

We derive square-root based algorithms for structured matrices and discuss potential applications to interpolation and matrix completion problems. The mathematical machinery used here is based on a standard Gaussian elimination technique and on simple results from matrix and linear system theory. We show how to exploit the inherent displacement structure in order to construct a convenient transmission-line cascade that makes evident the required interpolation and completion conditions. We also introduce the concept of time-variant structured matrices and discuss its applications to matrix completions and time-variant interpolation problems.

This paper continues some recent work on the numerical solution of the steady incompressible Navier–Stokes equations. We present a new method, similar to the one presented in Rebholz et al., but with superior convergence and numerical properties. The method is efficient as it allows one to solve the same symmetric positive‐definite system for the pressure at each iteration, allowing for the simple preconditioning and the reuse of preconditioners. We also demonstrate how one can replace the Schur complement system with a diagonal matrix inversion while maintaining accuracy and convergence, at a small fraction of the numerical cost. Convergence is analyzed for Newton and Picard‐type algorithms, as well as for the Schur complement approximation.

As a kind of commutative and associative four‐dimensional algebra, commutative quaternion has better applications in color image and signal processing than quaternion. Matrix decomposition is of great concern in the theoretical study and numerical calculation of commutative quaternion. Two kinds of disadvantages of commutative quaternion make the decomposition of commutative quaternion matrix extremely difficult. On one hand, commutative quaternion is not a kind of complete four dimensional division algebra because of the zero divisors. On the other hand, computing the inverse of commutative quaternion is very complicated. In this paper, the semi‐tensor product (STP) of matrices will be used to overcome the above two kinds of shortcomings. And we will propose a real structure‐preserving algorithm based on STP of matrices for commutative quaternion LU decomposition, which makes full use of high‐level operations, relation of operations between commutative quaternion matrices and their ‐representation matrices. Numerical experiments will be provided to demonstrate the efficiency of the real structure‐preserving algorithm based on STP of matrices. Meanwhile, we will apply the structure‐preserving algorithm for strict image authentication.