Search

This paper addresses the mixed $\mathcal{H}_2$/$\mathcal{H}_\infty$ control problem for linear discrete-time systems with structured multiplicative input noises. The structured random noises are modelled by a diagonal matrix where individual multiplicative noise is allowed to be different. The Nash game methodology is adopted to deal with such a multiobjective control problem, and a necessary and sufficient condition is analytically given in terms of two cross-coupled modified algebraic Riccati equations (MAREs). Based on the stabilizing solutions to MAREs, bounded causal equilibrium strategies are thus conducted and the resulting control law optimizes some specific performance together with a guaranteed worst case performance. Some relevant issues that can be regarded as special cases of the problem under consideration are also discussed. Moreover, a numerical example is included to show the validity of the present results.

We develop two fast algorithms for Hessenberg reduction of a structured matrix $A = D + UV^H$ where $D$ is a real or unitary $n \times n$ diagonal matrix and $U, V \in\mathbb{C}^{n \times k}$. The proposed algorithm for the real case exploits a two--stage approach by first reducing the matrix to a generalized Hessenberg form and then completing the reduction by annihilation of the unwanted sub-diagonals. It is shown that the novel method requires $O(n^2k)$ arithmetic operations and it is significantly faster than other reduction algorithms for rank structured matrices. The method is then extended to the unitary plus low rank case by using a block analogue of the CMV form of unitary matrices. It is shown that a block Lanczos-type procedure for the block tridiagonalization of $\Re(D)$ induces a structured reduction on $A$ in a block staircase CMV--type shape. Then, we present a numerically stable method for performing this reduction using unitary transformations and we show how to generalize the sub-diagonal elimination to this shape, while still being able to provide a condensed representation for the reduced matrix. In this way the complexity still remains linear in $k$ and, moreover, the resulting algorithm can be adapted to deal efficiently with block companion matrices., 25 pages

We present a new tunably accurate Laguerre Petrov--Galerkin spectral method for solving linear multiterm fractional initial value problems with derivative orders at most one and constant coefficients on the half line. Our method results in a matrix equation of special structure which can be solved in $\mathcal{O}(N \log N)$ operations. We also take advantage of recurrence relations for the generalized associated Laguerre functions (GALFs) in order to derive explicit expressions for the entries of the stiffness and mass matrices, which can be factored into the product of a diagonal matrix and a lower-triangular Toeplitz matrix. The resulting spectral method is efficient for solving multiterm fractional differential equations with arbitrarily many terms, which we demonstrate by solving a fifty-term example. We apply this method to a distributed order differential equation, which is approximated by linear multiterm equations through the Gauss--Legendre quadrature rule. We provide numerical examples demonstra...

Laplacian pyramid based Laurent polynomial (LP$^2$) matrices are generated by Laurent polynomial column vectors and have long been studied in connection with Laplacian pyramidal algorithms in Signal Processing. In this paper, we investigate when such matrices are scalable, that is when right multiplication by Laurent polynomial diagonal matrices results in paraunitary matrices. The notion of scalability has recently been introduced in the context of finite frame theory and can be considered as a preconditioning method for frames. This paper significantly extends the current research on scalable frames to the setting of polyphase representations of filter banks. Furthermore, as applications of our main results we propose new construction methods for tight wavelet filter banks and tight wavelet frames., Version accepted for publication in SIAM Journal on Matrix Analysis and Applications

We prove a new inequality for the expectation $E\left[\log\det\left(\mathbf{W}\mathbf{Q}+\mathbf{I}\right)\right]$, where $\mathbf{Q}$ is a nonnegative definite matrix and $\mathbf{W}$ is a diagonal random matrix with identically distributed nonnegative diagonal entries. A sharp lower bound is obtained by substituting $\mathbf{Q}$ by the diagonal matrix of its eigenvalues $\mathbf{\Gamma}$. Conversely, if this inequality holds for all $\mathbf{Q}$ and $\mathbf{\Gamma}$, then the diagonal entries of $\mathbf{W}$ are necessarily identically distributed. From this general result, we derive related deterministic inequalities of Muirhead- and Rado-type. We also present some applications in information theory: We derive bounds on the capacity of parallel Gaussian fading channels with colored additive noise and bounds on the achievable rate of noncoherent Gaussian fading channels.

This work is motivated by two two-sided optimization problems that arise in atomic chemistry when one is looking for a minimal set of localized orbitals with prescribed occupation numbers, respectively, a prescribed total number of electrons. We first propose an optimal analytic solution of the first problem, and then show that an optimal solution of the second problem can be found by solving a convex quadratic programming problem with box constraints and $p$ unknowns. We prove that the latter problem is solved by the active-set method in at most $2p$ iterations. These investigations reveal that the optimal solutions of both problems are projections that are in $\mathcal{C}=\{Y\in\mathbb{R}^{n\times p}:Y^TY=I_p, Y^TN Y=\Delta\}$ for $N$ symmetric and $\Delta$ diagonal and that these solutions are generally nonunique. The question of how one can then select a particular solution out of the set $\mathcal{C}$ leads to the main problem of this paper: the optimization of an arbitrary smooth function whose mini...

Limited-memory incomplete Cholesky factorizations can provide robust preconditioners for sparse symmetric positive-definite linear systems. In this paper, the focus is on extending the approach to sparse symmetric indefinite systems in saddle-point form. A limited-memory signed incomplete Cholesky factorization of the form $LDL^T$ is proposed, where the diagonal matrix $D$ has entries $\pm 1$. The main advantage of this approach is its simplicity as it avoids the use of numerical pivoting. Instead, a global shift strategy involving two shifts (one for the $(1,1)$ block and one for the (2,2) block of the saddle-point matrix) is used to prevent breakdown and to improve performance. The matrix is optionally prescaled and preordered using a standard sparse matrix ordering scheme that is then postprocessed to give a constrained ordering that reduces the likelihood of breakdown and need for shifts. The use of intermediate memory (memory used in the construction of the incomplete factorization but subsequently d...

We derive a new representation of Lagrangian subspaces in the form ${\rm Im} \Pi^T\big[\begin{smallmatrix}I \\ X\end{smallmatrix}\big]$, where $\Pi$ is a symplectic matrix which is the product of a permutation matrix and a real orthogonal diagonal matrix, and $X$ satisfies $\left\vert X_{ij}\right\vert \leq \begin{cases}1 & \text{if $i=j$,}\\ \sqrt{2} & \text{if $i\neq j$.} \end{cases}$ This representation allows us to limit element growth in the context of doubling algorithms for the computation of Lagrangian subspaces and the solution of Riccati equations. It is shown that a simple doubling algorithm using this representation can reach full machine accuracy on a wide range of problems, obtaining invariant subspaces of the same quality as those computed by the state-of-the-art algorithms based on orthogonal transformations. The same idea carries over to representations of arbitrary subspaces and can be used for other types of structured pencils.

This paper presents two methods based on domain decomposition concepts for determining the diagonal of the inverse of specific matrices. The first uses a divide-and-conquer principle and the Sherman-Morrison-Woodbury formula and assumes that the matrix can be decomposed into a $2 \times 2$ block-diagonal matrix and a low-rank matrix. The second method is a standard domain decomposition approach in which local solves are combined with a global correction. Both methods can be successfully combined with iterative solvers and sparse approximation techniques. The efficiency of the methods usually depends on the specific implementation, which should be fine-tuned for different test problems. Preliminary results for some two-dimensional (2D) problems are reported to illustrate the proposed methods.

We consider the solutions of Hermitian positive definite Toeplitz-plus-diagonal systems $(T+D)x=b$, where $T$ is a Toeplitz matrix and $D$ is diagonal and positive. However, unlike the case of Toeplitz systems, no fast direct solvers have been developed for solving them. In this paper, we employ the preconditioned conjugate gradient method with approximate inverse circulant-plus-diagonal preconditioners to solving such systems. The proposed preconditioner can be constructed and implemented efficiently using fast Fourier transforms. We show that if the entries of $T$ decay away exponentially from the main diagonals, the preconditioned conjugate gradient method applied to the preconditioned system converges very quickly. Numerical examples including spatial regularization for image deconvolution application are given to illustrate the effectiveness of the proposed preconditioner.

We characterize ordered orthonormal bases $\{s_1,\dots,s_n\} \subset \mathbb{R}^n$ for which the matrix $s_1s_1^T+\dots+s_ts_t^T$ is nonnegative for each $t=1,\dots,n$. In particular we show that the associated orthogonal matrix $S=[s_1,s_2,\dots,s_n]$ has certain submatrices that are Soules matrices. Based on the construction of Soules bases developed by Elsner, Nabben, and Neumann [Linear Algebra Appl., 271 (1998), pp. 323-343] we are able to construct all possible such bases. Additionally, we show that this set is the closure of the set of Soules bases. We also include some results regarding matrix functions of $S \Lambda S^T$, where $\Lambda$ is a diagonal matrix with nonincreasing diagonal entries.

For a symmetric positive semidefinite diagonally dominant matrix, if its off-diagonal entries and its diagonally dominant parts for all rows (which are defined for a row as the diagonal entry subtracted by the sum of absolute values of off-diagonal entries in that row) are known to a certain relative accuracy, we show that its eigenvalues are known to the same relative accuracy. Specifically, we prove that if such a matrix is perturbed in a way that each off-diagonal entry and each diagonally dominant part have relative errors bounded by some $\epsilon$, then all its eigenvalues have relative errors bounded by $\epsilon$. The result is extended to the generalized eigenvalue problem.

In this paper we introduce a new orthogonalization method. Given a real $m \times n$ matrix $A$, the new method constructs an SVD-type decomposition of the form $A = \hat U\hat\Sigma\hat V^T$. The columns of $\hat U$ and $\hat V$ are orthonormal, or nearly orthonormal, while $\hat\Sigma$ is a diagonal matrix whose diagonal entries approximate the singular values of $A$. The method has three versions: a “left-side" orthogonalization scheme in which the columns of $\hat U$ constitute an orthonormal basis of Range$(A)$, a “right-side" orthogonalization scheme in which the columns of $\hat V$ constitute an orthonormal basis of Range$(A^T)$, and a third version in which both $\hat U$ and $\hat V$ have orthonormal columns, but the decomposition is not exact. The new decompositions may substitute the SVD in many applications.

An $n \times n$ real matrix $A$ is said to be (totally negative) totally nonpositive if every minor is (negative) nonpositive. In this paper, we study the properties of a totally nonpositive matrix and characterize the case of a nonsingular totally nonpositive matrix $A$, with $a_{11}< 0$ in terms of its $LDU$ factorization ($L(U)$) is a unit lower- (upper-) triangular matrix, respectively, and $D$ is a diagonal matrix). This characterization allows us to significantly reduce the number of minors to be checked in order to decide the total nonpositivity of a nonsingular matrix with a negative (1,1) entry.

We derive absolute perturbation bounds for the coefficients of the characteristic polynomial of a $n\times n$ complex matrix. The bounds consist of elementary symmetric functions of singular values, and suggest that coefficients of normal matrices are better conditioned with regard to absolute perturbations than those of general matrices. When the matrix is Hermitian positive-definite, the bounds can be expressed in terms of the coefficients themselves. We also improve absolute and relative perturbation bounds for determinants. The basis for all bounds is an expansion of the determinant of a perturbed diagonal matrix.

Let $(X_n)_{n\ge 1}$ be a sequence of independent centered random vectors in ${\bf R}^d$. We give conditions under which the sequence $S_n=\sum_{i=1}^nX_i$ normalized by a matricial sequence $(H_n)$ satisfies a compact law of the iterated logarithm. As an application of this result, we obtain the compact law of the iterated logarithm for $B_n^{-1/2}S_n$ and for $\Delta_n^{-1/2}S_n$, where $B_n$ is the covariance matrix of $S_n$, and where $\Delta_n$ is the diagonal matrix whose jth diagonal term is the jth diagonal term of $B_n$; the eigenvalues of $B_n$ may go to infinity with different rates, but their iterated logarithms have to be equivalent.

Explicit expressions for the entries of the inverse of a general diagonal matrix are derived via the decomposition approach. In particular, three-term recurrence relations with variable coefficients are given, and with their solutions closed form expressions for each entry of the fundamental matrix of a Markov chain-random walk (MC-RW) model are established. We also obtain explicit formulas for the absorption probabilities and the mean time to absorption as well as the mean number of steps before absorption of the MC-RW model in the presence of two semiabsorbing boundaries.

In this paper we present a fast direct solver for certain classes of dense structured linear systems that works by first converting the given dense system to a larger system of block sparse equations and then uses standard sparse direct solvers. The kind of matrix structures that we consider are induced by numerical low rank in the off-diagonal blocks of the matrix and are related to the structures exploited by the fast multipole method (FMM) of Greengard and Rokhlin. The special structure that we exploit in this paper is captured by what we term the hierarchically semiseparable (HSS) representation of a matrix. Numerical experiments indicate that the method is probably backward stable.

Sufficient conditions for global convergence of the primal-dual active set strategy for finite and infinite dimensional quadratic, as well as nonlinear optimization, problems with affine equality and inequality constraints are presented. These conditions involve diagonal dominance and cone preserving properties of the operator defining the cost functional. Globalization strategies are also provided, and specific sufficient conditions for the primal-dual active set step to have a descent property are given.

In this paper we study the existence of linear quadratic Gaussian (LQG)-balanced realizations for continuous-time infinite-dimensional systems. LQG-balanced realizations are those for which the optimal cost operator for the system and its dual system are equal (and diagonal). The class of systems we consider is that of distributional resolvent linear systems which includes well-posed linear systems as a subclass. We prove the existence of LQG-balanced realizations under a finite cost condition for both the system and its dual system. We also show that an LQG-balanced realization of a well-posed transfer function is well-posed. We further show that approximately controllable and observable LQG-balanced realizations are unique up to a unitary state-space transformation. Finally, we show that the spectrum of the product of the optimal cost operator of a system and its dual system is independent of the particular realization. Our method of proof shows the connections with coprime factorizations, Lyapunov-balanced realizations, and discrete-time systems. The main reason for studying LQG-balanced realizations is that truncated LQG-balanced realizations provide a good approximation of the original system. We show that, under certain conditions, this is also true in the infinite-dimensional case by proving an error bound in the gap-metric.

We consider the problem of finding a common quadratic Lyapunov function to demonstrate the stability of a not necessarily bounded family of matrices which incorporate design freedoms. Generically, this can be viewed as the problem of picking a family of controller (or observer) gains so that the family of closed-loop system matrices admits a common quadratic Lyapunov function. We provide several necessary and sufficient conditions for various structures of matrix families. Specifically, we consider matrix families which can be obtained through similarity transformations from a lower Hessenberg structure. Families of matrices containing a subset of diagonal matrices, which is invariant under the design freedoms, are also considered since they occur in many applications. The conditions for uniform solvability of the Lyapunov inequalities are explicitly given and involve inequalities regarding relative magnitudes of terms in the matrices. Various motivating applications of the obtained results to observer and controller designs for time-varying, switched, and nonlinear systems are highlighted.

The four existing stable factorization methods for symmetric indefinite matrices suffer serious defects when applied to banded matrices. Partial pivoting (row or column exchanges) maintains a band structure in the reduced matrix and the factors, but destroys symmetry completely once an off-diagonal pivot is used. Two-by-two block pivoting and Gaussian reduction to tridiagonal (Aasen's algorithm) maintain symmetry at all times, but quickly destroy the band structure in the reduced matrices. Orthogonal reductions to tridiagonal maintain both symmetry and the band structure, but are too expensive for linear-equation solvers. We propose a new pivoting method, which we call snap-back pivoting. When applied to banded symmetric matrices, it maintains the band structure (like partial pivoting does), it keeps the reduced matrix symmetric (like 2-by-2 pivoting and reductions to tridiagonal), and it is fast. Snap-back pivoting reduces the matrix to a diagonal form using a sequence of elementary elimination steps, most of which are applied symmetrically from the left and from the right (but some are applied unsymmetrically). In snap-back pivoting, if the next diagonal element is too small, the next pivoting step might be unsymmetric, leading to asymmetry in the next row and column of the factors. But the reduced matrix snaps back to symmetry once the next step is completed.

A matrix $A\in M_{n}(\mathbb{C})$ is called unitoid if it is congruent to a diagonal matrix. Necessary and sufficient conditions are given on the canonical angles of a unitoid matrix so that sufficiently small perturbations remain unitoid. This, in particular, resolves the question of when simultaneous diagonalizability of two Hermitian matrices is retained under perturbation.

In this paper we investigate some matrix structures on $\cee^{n\times n}$ that have a good behavior under matrix inversion. The first type of structure is closely related to low displacement rank matrices. Next, we show that for a matrix having a low rank submatrix, the inverse matrix also must have a low rank submatrix, which we can explicitly determine. This allows us to generalize a theorem due to Fiedler and Markham. The generalization consists in the fact that our rank structures may have a certain correction term, which we call the shift matrix $\Lam_k\in\mathbb{C}^{m \times m}$, for suitable m, and with Fiedler and Markham's theorem corresponding to the limiting cases $\Lam_k\to 0$ and $\Lam_k\to\infty I$.

We present a new diagonal balancing technique for regular matrix pencils $\lambda B-A$, which aims at reducing the sensitivity of the corresponding generalized eigenvalues. It is inspired by the balancing technique of a square matrix A and has a comparable complexity. The diagonally scaled pencil has row and column norms that are balanced in a precise sense. We also show that balancing a pencil boils down to making it closer to some standardized normal pencil. We give numerical examples illustrating that the sensitivity of generalized eigenvalues of a pencil may significantly improve after balancing.

In this paper we investigate some matrix structures on $\mathbb{C}^{m\times n}$ that are preserved by Schur complementation. The first type of structure is closely related to low displacement rank matrices. Next, we show that for a matrix having a low rank submatrix, the Schur complement also must have a low rank submatrix, which we can explicitly determine. This property holds even if the low rank submatrix contains a certain correction term that we call the shift matrix.

In this paper we study the variation of the spectrum of block-diagonal systems under perturbations of compatible block structure with fixed zero blocks at arbitrarily prescribed locations ("Gershgorin-type perturbations"). We derive explicit and computable formulae for the associated mu-values. The results are then applied to characterize spectral value sets and stability radii for such perturbed systems. By specializing our results to the scalar diagonal case, the classical eigenvalue inclusion theorems of Gershgorin, Brauer, and Brualdi are obtained as corollaries. Moreover it follows that the inclusion regions of Brauer and Brualdi are optimal for the corresponding perturbation structures.

Consider solving a linear program in standard form where the constraint matrix $A$ is $m \times n$, with $n \gg m \gg 1$. Such problems arise, for example, as the result of finely discretizing a semi-infinite program. The cost per iteration of typical primal-dual interior-point methods on such problems is $O(m^2n)$. We propose to reduce that cost by replacing the normal equation matrix, $AD^2A^{\T}$, where $D$ is a diagonal matrix, with a "reduced" version (of same dimension), $A_QD_Q^2A_Q^{\T}$, where $Q$ is an index set including the indices of $M$ most nearly active (or most violated) dual constraints at the current iterate, with $M\geq m$ a prescribed integer. This can result in a speedup of close to $n/|Q|$ at each iteration. Promising numerical results are reported for constraint-reduced versions of a dual-feasible affine-scaling algorithm and of Mehrotra's predictor-corrector method [S. Mehrotra, it SIAM J. Optim., 2 (1992), pp. 575-601]. In particular, while it could be expected that neglecting a large portion of the constraints, especially at early iterations, may result in a significant deterioration of the search direction, it appears that the total number of iterations typically remains essentially constant as the size of the reduced constraint set is decreased down to some threshold. In some cases this threshold is a small fraction of the total set. In the case of the affine-scaling algorithm, global convergence and local quadratic convergence are proved.

An extended sequentially semiseparable (SSS) representation derived from time-varying system theory is used to capture, on the one hand, the low-rank of the off-diagonal blocks of a matrix for the purposes of efficient computations and, on the other, to provide for sufficient descriptive richness to allow for backward stability in the computations. We present (i) a fast algorithm (linear in the number of equations) to solve least squares problems in which the coefficient matrix is in SSS form, (ii) a fast algorithm to find the SSS form of X such that AX=B, where A and B are in SSS form, and (iii) a fast model reduction technique to improve the SSS form.

For an $n \times n$ tridiagonal matrix we exploit the structure of its QR factorization to devise two new algorithms for computing the 1-norm condition number in $O(n)$ operations. The algorithms avoid underflow and overflow, and are simpler than existing algorithms since tests are not required for degenerate cases. An error analysis of the first algorithm is given, while the second algorithm is shown to be competitive in speed with existing algorithms. We then turn our attention to an $n \times n$ diagonal-plus-semiseparable matrix, $A$, for which several algorithms have recently been developed to solve $Ax=b$ in $O(n)$ operations. We again exploit the QR factorization of the matrix to present an algorithm that computes the 1-norm condition number in $O(n)$ operations.

Let $B$ be an $m\times n$ ($m\ge n$) complex (or real) matrix. It is known that there is a unique {\em polar decomposition} $B=QH$, where $Q^*Q=I$, the $n\times n$ identity matrix, and $H$ is positive definite, provided that $B$ has full column rank. If $B$ is perturbed to $\wtd B$, how do the polar factors $Q$ and $H$ change? This question has been investigated quite extensively, but most work so far has been on how the perturbation changed the unitary polar factor $Q$, with very little on the positive polar factor $H$, except $\|H-\wtd H\|_{\F}\le\sqrt 2\|B-\wtd B\|_{\F}$ in the Frobenius norm, due to [F. Kittaneh, {\em Comm.\ Math.\ Phys.}, 104 (1986), pp. 307--310], where $\wtd Q$ and $\wtd H$ are the corresponding polar factors of $\wtd B$. While this inequality of Kittaneh shows that $H$ is always well behaved under perturbations, it does not tell much about smaller entries of $H$ in the case when $H$'s entries vary a great deal in magnitude. This paper is intended to fill the gap by addressing the variations of $H$ for a graded matrix $B=GS$, where $S$ is a scaling matrix and usually diagonal (but may not be). The elements of $S$ can vary wildly, while $G$ is well conditioned. In such cases, the magnitudes of $H$'s entries indeed often vary a lot, and thus any bound on $\|H-\wtd H\|_{\F}$ means little, if anything, to the accuracy of $\wtd H$'s smaller entries. This paper proposes a new way of measuring the errors in the $H$ factor via bounding the scaled difference $(\wtd H-H)S^{-1}$, as well as accurately computing the factor when $S$ {\em is\/} diagonal. Numerical examples are presented. The results are also extended to the matrix square root of a graded positive definite matrix.

In this paper, we study factorized banded inverse preconditioners for matrices with Toeplitz structure. We show that if a Toeplitz matrix T has certain off-diagonal decay property, then the factorized banded inverse preconditioner approximates T-1 accurately, and the spectra of these preconditioned matrices are clustered around 1. In nonlinear image restoration applications, Toeplitz-related systems of the form I + T* D T are required to solve, where D is a positive nonconstant diagonal matrix. We construct inverse preconditioners for such matrices. Numerical results show that the performance of our proposed preconditioners are superior to that of circulant preconditioners. A two-dimensional nonlinear image restoration example is also presented to demonstrate the effectiveness of the proposed preconditioner.

This paper presents a preconditioning method based on combining two-sided permutations with a multilevel approach. The nonsymmetric permutation exploits a greedy strategy to put large entries of the matrix in the diagonal of the upper leading submatrix. The method can be regarded as a complete pivoting version of the incomplete LU factorization. This leads to an effective incomplete factorization preconditioner for general nonsymmetric, irregularly structured, sparse linear systems.

We consider the Gersgorin disc separation from the origin for (doubly) diagonally dominant matrices and their Schur complements, showing that the separation of the Schur complement of a (doubly) diagonally dominant matrix is greater than that of the original grand matrix. As application we discuss the localization of eigenvalues and present some upper and lower bounds for the determinant of diagonally dominant matrices.

In this paper we present a new iteration-complexity bound for the Mizuno--Todd--Ye predictor-corrector (MTY P-C) primal-dual interior-point algorithm for linear programming. The analysis of the paper is based on the important notion of crossover events introduced by Vavasis and Ye. For a standard form linear program min{cTx:Ax = b, x \ge 0} with decision variable $x\in \Re^n$, e show that the MTY P-C algorithm, started from a well-centered interior-feasible solution with duality gap $n \mu_0$, finds an interior-feasible solution with duality gap less than $n\eta$ in ${\cal O}(T(\mu_0/\eta)+n^{3.5}\log({{\bar\chi^*_A}}))$ iterations, where $T(t)\equiv\min\{n^2\log(\log t) ,\log t\}$ for all $t>0$ and ${{\bar\chi^*_A}}$ is a scaling invariant condition number associated with the matrix $A$. More specifically, ${{\bar\chi^*_A}}$ is the infimum of all the conditions numbers ${\bar\chi}_{AD}$, where D varies over the set of positive diagonal matrices. Under the setting of the Turing machine model, our analysis yields an ${\cal O}(n^{3.5}L_A+ \min\{n^2\log L,L\})$ iteration-complexity bound for the MTY P-C algorithm to find a primal-dual optimal solution, where LA and L are the input sizes of the matrix A and the data (A,b,c), respectively. This contrasts well with the classical iteration-complexity bound for the MTY P-C algorithm, which depends linearly on L instead of log L.

Solving systems of linear equations with "normal" matrices of the form A D2 AT is a key ingredient in the computation of search directions for interior-point algorithms. In this article, we establish that a well-known basis preconditioner for such systems of linear equations produces scaled matrices with uniformly bounded condition numbers as D varies over the set of all positive diagonal matrices. In particular, we show that when A is the node--arc incidence matrix of a connected directed graph with one of its rows deleted, then the condition number of the corresponding preconditioned normal matrix is bounded above by m(n - m + 1), where m and n are the number of nodes and arcs of the network.

This paper is concerned with the characterization of the Hurwitz stability of matrices, which can be written as the product of two square matrices AD with A precisely known and D belonging to the set of all positive diagonal matrices, and the Schur stability of matrices AD for all diagonal D, whose entries have absolute magnitude less than or equal to 1, known as the problem of ${\Bbb D}$-stability. Sufficient conditions are given in terms of linear matrix inequalities formulated at the vertices of an adequately chosen polytope domain, allowing simple and numerically efficient evaluations of ${\Bbb D}$-stability. The conditions proposed provide less conservative results and encompass previous conditions from the literature, as illustrated by examples.

The Davidson method is a popular technique to compute a few of the smallest (or largest) eigenvalues of a large sparse real symmetric matrix. It is effective when the matrix is nearly diagonal, that is, when the matrix of eigenvectors is close to the identity matrix. However, its convergence properties are not yet well understood, and neither is how it behaves compared to the more recent Jacobi--Davidson method, for which a proper convergence analysis exists. In this paper, we develop a new convergence analysis of the Davidson method. This analysis proves that the convergence is fast for nearly diagonal matrices when the method is initialized in the standard way. One may at this stage not expect any significant improvement by shifting to the Jacobi--Davidson method. On the other hand, the latter may be more effective for more general initial approximations. It is also best suited for matrices that are not nearly diagonal, thanks to the use of more sophisticated preconditioning and/or inner iterations.

In principle, the exponential of a matrix could be computed in many ways. Methods involving approximation theory, differential equations, the matrix eigenvalues, and the matrix characteristic polyn...

Limited-memory BFGS quasi-Newton methods approximate the Hessian matrix of second derivatives by the sum of a diagonal matrix and a fixed number of rank-one matrices. These methods are particularly effective for large problems in which the approximate Hessian cannot be stored explicitly. It can be shown that the conventional BFGS method accumulates approximate curvature in a sequence of expanding subspaces. This allows an approximate Hessian to be represented using a smaller reduced matrix that increases in dimension at each iteration. When the number of variables is large, this feature may be used to define limited-memory reduced-Hessian methods in which the dimension of the reduced Hessian is limited to save storage. Limited-memory reduced-Hessian methods have the benefit of requiring half the storage of conventional limited-memory methods. In this paper, we propose a particular reduced-Hessian method with substantial computational advantages compared to previous reduced-Hessian methods. Numerical results from a set of unconstrained problems in the CUTE test collection indicate that our implementation is competitive with the limited-memory codes L-BFGS and L-BFGS-B.

An iterative method for identifying generalized diagonally dominant matrices\break (GDDMs, or H-matrices) was given in [B. Li et al., Linear Algebra Appl., 271 (1998), pp. 179--190], where the method is divergent when the matrix is not a GDDM. In this paper, we present an improved version. The new method is always convergent and needs fewer iterations than the earlier one. Some interesting features of the new method are presented. Spectral radii of nonnegative matrices with a constant diagonal entry also can be computed by our method.

We consider the rapid simulation of three-dimensional electromagnetic problems in geophysical parameter regimes, where the conductivity may vary significantly and the range of frequencies is moderate. Toward developing a multigrid preconditioner, we present a Fourier analysis based on a finite-volume discretization of a vector potential formulation of time-harmonic Maxwell's equations on a staggered grid in three dimensions. We prove grid-independent bounds on the eigenvalue and singular value ranges of the system obtained using a preconditioner based on exact inversion of the dominant diagonal blocks of the non-Hermitian coefficient matrix. This result implies that a preconditioner that uses single multigrid cycles to effect inversion of the diagonal blocks also yields a preconditioned system with an $\ell_2$-condition number bounded independent of the grid size. We then present numerical examples for more realistic situations involving large variations in conductivity (i.e.,\ jump discontinuities). Block-preconditioning with one multigrid cycle using Dendy's BOXMG solver is found to yield convergence in very few iterations, apparently independent of the grid size. The experiments show that the somewhat restrictive assumptions of the Fourier analysis do not prohibit it from describing the essential local behavior of the preconditioned operator under consideration. A very efficient, practical solver is obtained.

Computing the linear least-squares estimate of a high-dimensional random quantity given noisy data requires solving a large system of linear equations. In many situations, one can solve this system efficiently using a Krylov subspace method, such as the conjugate gradient (CG) algorithm. Computing the estimation error variances is a more intricate task. It is difficult because the error variances are the diagonal elements of a matrix expression involving the inverse of a given matrix. This paper presents a method for using the conjugate search directions generated by the CG algorithm to obtain a convergent approximation to the estimation error variances. The algorithm for computing the error variances falls out naturally from a new estimation-theoretic interpretation of the CG algorithm. This paper discusses this interpretation and convergence issues and presents numerical examples. The examples include a 105-dimensional estimation problem from oceanography.

Sparse approximate inverses are considered as smoothers for multigrid. They are based on the SPAI-Algorithm [M. J. Grote and T. Huckle, SIAM J. Sci. Comput., 18 (1997), pp. 838--853], which constructs a sparse approximate inverse M of a matrix A by minimizing I -MA in the Frobenius norm. This yields a new hierarchy of smoothers: SPAI-0, SPAI-1, SPAI$(\varepsilon)$. Advantages of SPAI smoothers over classical smoothers are inherent parallelism, possible local adaptivity, and improved robustness. The simplest smoother, SPAI-0, is based on a diagonal matrix M. It is shown to satisfy the smoothing property for symmetric positive definite problems. Numerical experiments show that SPAI-0 smoothing is usually preferable to damped Jacobi smoothing. For the SPAI-1 smoother the sparsity pattern of M is that of A; its performance is typically comparable to that of Gauss--Seidel smoothing; however, both the computation and the application of the smoother remain inherently parallel. In more difficult situations, where the simpler SPAI-0 and SPAI-1 smoothers are not adequate, the SPAI$(\varepsilon)$ smoother provides a natural procedure for improvement where needed. Numerical examples illustrate the usefulness of SPAI smoothing.

We consider bipartite matching algorithms for computing permutations of a sparse matrix so that the diagonal of the permuted matrix has entries of large absolute value. We discuss various strategies for this and consider their implementation as computer codes. We also consider scaling techniques to further increase the relative values of the diagonal entries. Numerical experiments show the effect of the reorderings and the scaling on the solution of sparse equations by a direct method and by preconditioned iterative techniques.

Let $D$ be a diagonal matrix and $E_{ij}$ denote the $n$-by-$n$ matrix with a $1$ in entry $(i,j)$ and $0$ in every other entry. An $n$-by-$n$ matrix $A$ has a successively ordered elementary bidiagonal $(SEB)$ factorization if it can be factored as \begin{equation*} A= \left(\prod_{k=1}^{n-1} \prod_{j=n}^{k+1} L_j(s_{jk})\right) \; D \; \left(\prod_{k=n-1}^{1} \prod_{j=k+1}^{n}U_j(t_{kj})\right), \end{equation*} in which $L_j(s_{jk})=I+s_{jk}E_{j,j-1}$ and $U_j(t_{kj})=I+t_{kj}E_{j-1,j}$ for some scalars $s_{jk},t_{kj}$. Note that some of the parameters $% s_{jk},t_{kj}$ may be zero, and the order of the bidiagonal factors is fixed. If this factorization corresponds to reduction of $A$ to $D$ via successive row/column operations in the specified order, it is called an elimination $SEB$ factorization. New rank conditions are formulated that are proved to be necessary and sufficient for matrix $A$ to have such a factorization. These conditions are related to known but more restrictive properties that ensure a bidiagonal factorization as above, but with all parameters $s_{jk},t_{kj}$ nonzero.

In this paper, we present a comparison of four preconditioning strategies for Jacobian systems arising in the fully implicit solution of radiation diffusion coupled with material energy transfer. The four preconditioning methods are block Jacobi, Schur complement, and operator splitting approaches that split the preconditioner solve into two steps. One splitting method includes the coupling of the radiation and material fields that appears in the matrix diagonal in the first solve, and the other method puts this coupling into the second solve. All preconditioning approaches use multigrid methods to invert blocks of the matrix formed from the diffusion operator. The Schur complement approach is clearly seen to be the most effective for a large range of weightings between the diffusion and energy coupling terms. In addition, tabulated opacity studies were conducted where, again, the Schur preconditioner performed well. Last, a parallel scaling study was done showing algorithmic scalability of the Schur preconditioner.

This paper presents a new algorithm for accurate floating-point computation of the singular value decomposition (SVD) of the product $A = B^{\tau} S C,$ where $ B\in {\hbox{{\bf R}$^{p\times m} $}},$ $C\in {\hbox{{\bf R}$^{q\times n}$}},$ $S\in {\hbox{{\bf R}$^{p\times q}$}},$ and $p\leq m,$ $q\leq n$.\ The new algorithm uses diagonal scalings, the LU factorization with complete pivoting, the QR factorization with column pivoting, and matrix multiplication to replace $A$ by $A' = B'^{\tau}S'C',$ where A and A' have the same singular values and the matrix A' is computed explicitly. The singular values of A' are computed using the Jacobi SVD algorithm. It is shown that the accuracy of the new algorithm is determined by (i) the accuracy of the QR factorizations of $B^{\tau}$ and $C^{\tau}$; (ii) the accuracy of the LU factorization with complete pivoting of S; and (iii) the accuracy of the computation of the SVD of a matrix $A'$ with moderate $\min_{D=\diag} \kappa_2(A'D)$.\ Theoretical analysis and numerical evidence show that, in the case of rank (B)= rank(C)=p and full rank S, the accuracy of the new algorithm is unaffected by replacing B, S, C with, respectively, D1 B, D2SD3,D4 C, where Di, i=1, . . .,4, are arbitrary diagonal matrices. As an application, the paper proposes new accurate algorithms for computing the (H,K)--SVD and (H-1,K)--SVD of S.

Special ultrametric matrices are, in a sense, extremal matrices in the boundary of the set of ultrametric matrices introduced by Martinez, Michon, and San Martin [ SIAM J. Matrix Anal. Appl., 15 (1994), pp. 98--106]. We show a simple construction of these matrices, if of order n, from nonnegatively edge-weighted trees on n vertices, or, equivalently, from nonnegatively edge-weighted paths. A general ultrametric matrix is then the sum of a nonnegative diagonal matrix and a special ultrametric matrix, with certain conditions fulfilled. The rank of a special ultrametric matrix is also recognized and it is shown that its Moore--Penrose inverse is a generalized diagonally dominant M-matrix. Some results on the nonsymmetric case are included.

This paper deals with convergence estimates of GMRES(m) [Saad and Schultz, { SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856--869] preconditioned by multigrid [Brandt, Math. Comp., 31 (1977), pp. 333--390], [Hackbusch, Multi-Grid Methods and Applications, Springer, Berlin, 1985]. Fourier analysis is a well-known and useful tool in the multigrid community for the prediction of two-grid convergence rates [Brandt, Math. Comp., 31 (1977), pp. 333--390], [Stuben and Trottenberg, in Multigrid Methods, Lecture Notes in Math. 960, K. Stuben and U. Trottenberg, eds., Springer, Berlin, pp. 1--176]. This analysis is generalized here to the situation in which multigrid is a preconditioner, since it is possible to obtain the whole spectrum of the two-grid iteration matrix. A preconditioned Krylov subspace acceleration method like GMRES(m) implicitly builds up a minimal residual polynomial. The determination of the polynomial coefficients is easily possible and can be done explicitly since, from Fourier analysis, a simple block-diagonal two-grid iteration matrix results. Based on the GMRES(m) polynomial, sharp theoretical convergence estimates can be obtained which are compared with estimates based on the spectrum of the iteration matrix. Several numerical scalar test problems are computed in order to validate the theoretical predictions.