Perturbation analysis of the orthogonal procrustes problem

Given two arbitrary real matricesA andB of the same size, the orthogonal Procrustes problem is to find an orthogonal matrixM such that the Frobenius norm MA – B is minimized. This paper treats the common case when the orthogonal matrixM is required to have a positive determinant. The stability of the problem is studied and supremum results for the perturbation bounds are derived.     

On condition numbers and the distance to the nearest ill-posed problem

Summary The condition number of a problem measures the sensitivity of the answer to small changes in the input. We call the problem ill-posed if its condition number is infinite. It turns out that for many problems of numerical analysis, there is a simple relationship between the condition number of a problem and the shortest distance from that problem to an ill-posed one: the shortest distance is proportional to the reciprocal of the condition number (or bounded by the reciprocal of the condition number). This is true for matrix inversion, computing eigenvalues and eigenvectors, finding zeros of polynomials, and pole assignment in linear control systems. In this paper we explain this phenomenon by showing that in all these cases, the condition number satisfies one or both of the diffrential inequalitiesm·²DM·², where D is the norm of the gradient of . The lower bound on D leads to an upper bound 1/m(x) on the distance. fromx to the nearest ill-posed problem, and the upper bound on D leads to a lower bound 1/(M(X)) on the distance. The attraction of this approach is that it uses local information (the gradient of a condition number) to answer a global question: how far away is the nearest ill-posed problem? The above differential inequalities also have a simple interpretation: they imply that computing the condition number of a problem is approximately as hard as computing the solution of the problem itself. In addition to deriving many of the best known bounds for matrix inversion, eigendecompositions and polynomial zero finding, we derive new bounds on the distance to the nearest polynomial with multiple zeros and a new perturbation result on pole assignment.     

Component-wise perturbation analysis and error bounds for linear least squares solutions

Perturbation bounds for the linear least squares problem min _x Ax –b₂ corresponding tocomponent-wise perturbations in the data are derived. These bounds can be computed using a method of Hager and are often much better than the bounds derived from the standard perturbation analysis. In particular this is true for problems where the rows ofA are of widely different magnitudes. Generalizing a result by Oettli and Prager, we can use the bounds to compute a posteriori error bounds for computed least squares solutions.     

Approximation of eigenfunctions in kernel-based spaces

Kernel-based methods in Numerical Analysis have the advantage of yielding optimal recovery processes in the “native” Hilbert space \(\mathcal {H}\) in which they are reproducing. Continuous kernels on compact domains have an expansion into eigenfunctions that are both L ₂-orthonormal and orthogonal in \(\mathcal {H}\) (Mercer expansion). This paper examines the corresponding eigenspaces and proves that they have optimality properties among all other subspaces of \(\mathcal {H}\). These results have strong connections to n-widths in Approximation Theory, and they establish that errors of optimal approximations are closely related to the decay of the eigenvalues. Though the eigenspaces and eigenvalues are not readily available, they can be well approximated using the standard n-dimensional subspaces spanned by translates of the kernel with respect to n nodes or centers. We give error bounds for the numerical approximation of the eigensystem via such subspaces. A series of examples shows that our numerical technique via a greedy point selection strategy allows to calculate the eigensystems with good accuracy.     

Generalizations of the projection method with applications to SOR theory for hermitian positive semidefinite linear systems

Summary Ann×n complex matrixB is calledparacontracting if B₂1 and 0x[N(I-B)]Bx₂<x₂. We show that a productB=B _k B _k–1 ...B ₁ ofk paracontracting matrices is semiconvergent and give upper bounds on the subdominant eigenvalue ofB in terms of the subdominant singular values of theB _i 's and in terms of the angles between certain subspaces. Our results here extend earlier results due to Halperin and due to Smith, Solomon and Wagner. We also determine necessary and sufficient conditions forn numbers in the interval [0, 1] to form the spectrum of a product of two orthogonal projections and hence characterize the subdominant eigenvalue of such a product. In the final part of the paper we apply the upper bounds mentioned earlier to provide an estimate on the subdominant eigenvalue of the SOR iteration matrix associated with ann×n hermitian positive semidefinite matrixA none of whose diagonal entries vanish.The work of this author was supported in part by NSF Research Grant No. MCS-8400879     

A superfast algorithm for multi-dimensional Padé systems

For a vector ofk+1 matrix power series, a superfast algorithm is given for the computation of multi-dimensional Padé systems. The algorithm provides a method for obtaining matrix Padé, matrix Hermite Padé and matrix simultaneous Padé approximants. When the matrix power series is normal or perfect, the algorithm is shown to calculate multi-dimensional matrix Padé systems of type (n ₀,...,n _k ) inO(n · log²n) block-matrix operations, where n=n ₀+...+n _k . Whenk=1 and the power series is scalar, this is the same complexity as that of other superfast algorithms for computing Padé systems. Whenk>1, the fastest methods presently compute these matrix Padé approximants with a complexity ofO(n²). The algorithm succeeds also in the non-normal and non-perfect case, but with a possibility of an increase in the cost complexity.Supported in part by NSERC grant No. A8035.Partially supported by NSERC operating grant No. 6194.     

Upper bounds for theL p -norms of martingales

Summary Let (f _n ) be a martingale. We establish a relationship between exponential bounds for the probabilities of the typeP(|f _n |>·T(f _n )) and the size of the constantC _p appearing in the inequality f ^* _p C _p T ^*(f) _p , for some quasi-linear operators acting on martingales.This research was supported in part by NSF Grant, no. DMS-8902418On leave from Academy of Physical Education, Warsaw, Poland     

A block incomplete orthogonalization method for large nonsymmetric eigenproblems

The incomplete orthogonalization method (IOM) proposed by Saad for computing a few eigenpairs of large nonsymmetric matrices is generalized into a block incomplete orthogonalization method (BIOM). It is studied how the departure from symmetry A – A ^H affects the conditioning of the block basis vectors generated by BIOM, and some relationships are established between the approximate eigenpairs obtained by BIOM and Ritz pairs. It is proved that BIOM behaves much like generalized block Lanczos methods if the basis vectors of the block Krylov subspace generated by it are strongly linearly independent. However, it is shown that BIOM may generate a nearly linearly dependent basis for a general nonsymmetric matrix. Numerical experiments illustrate the convergence behavior of BIOM.This work was supported in part by the Graduiertenkolleg at the University of Bielefeld, Germany.     

Extreme values and the law of the iterated logarithm

Summary IfX takes values in a Banach spaceB and is in the domain of attraction of a Gaussian law onB, thenX satisfies the compact law of the iterated logarithm (LIL) with respect to a regular normalizing sequence { _n } iffX satisfies a certain integrability condition. The integrability condition is equivalent to the fact that the maximal term of the sample {X ₁, X ₂,..., X _n} does not dominate the partial sums {S _n}, and here we examine the precise influence of these maximal terms and its relation to the compactLIL. In particular, it is shown that if one deletes enough of the maximal terms there is always a compactLIL with non-trivial limit set.Supported in part by NSF Grant MCS-8219742Work done while visiting the University of Wisconsin, Madison, with partial support by NSF Grant MCS-8219742     

Proper deflating subspaces: properties,algorithms and applications

An algorithm for computing proper deflating subspaces with specified spectrum for an arbitrary matrix pencil is presented. The method uses refined algorithms for computing the generalized Schur form of a matrix pencil and enlightens the connection that exists between reducing and proper deflating subspaces. The proposed algorithm can be applied for computing the stabilizing solution of the generalized algebraic Riccati equation, a recently introduced concept which extends the usual algebraic Riccati equation.     

On the norms of polynomials of two adjoint projections and a shift

LetP be a projection (non-selfadjoint in general), andV a selfadjoint involution acting in a Hilbert spaceH. In this paper the polynomialsF(X, Y, Z) of three non-commuting variables are described such that the norms F(P, P ^*,V) depend only on P. A method of calculation of the norms F(P, P ^*,V) for such polynomials is given. For polynomialsF(P, P ^*) this problem was investigated in [KMF], [FKM].     

Numerical stability of the Chebyshev method for the solution of large linear systems

Summary This paper contains the rounding error analysis for the Chebyshev method for the solution of large linear systemsAx+g=0 whereA=A ^* is positive definite. We prove that the Chebyshev method in floating point arithmetic is numerically stable, which means that the computed sequence {x _k} approximates the solution such that x _k – is of order AA ^–1 where is the relative computer precision.We also point out that in general the Chebyshev method is not well-behaved, which means that the computed residualsr _k=Ax _k+g are of order A²A ^–1.This work was supported in part by the Office of Naval Research under Contract N0014-67-0314-0010, NR 044-422 and by the National Science Foundation under Grant GJ32111     

Round-off error analysis of iterations for large linear systems

Summary We deal with the rounding error analysis of successive approximation iterations for the solution of large linear systemsA _x =b. We prove that Jacobi, Richardson, Gauss-Seidel and SOR iterations arenumerically stable wheneverA=A ^*>0 andA has PropertyA. This means that the computed resultx _k approximates the exact solution with relative error of order A·A ^–1 where is the relative computer precision. However with the exception of Gauss-Seidel iteration the residual vector Ax _k –b is of order A² A ^–1 and hence the remaining three iterations arenot well-behaved.This work was partly done during the author's visit at Carnegie-Mellon University and it was supported in part by the Office of Naval Research under Contract N00014-76-C-0370; NR 044-422 and by the National Science Foundation under Grant MCS75-222-55     

A Faster and Simpler Recursive Algorithm for the LAPACK Routine DGELS

We present new algorithms for computing the linear least squares solution to overdetermined linear systems and the minimum norm solution to underdetermined linear systems. For both problems, we consider the standard formulation min AX – B _F and the transposed formulation min A ^T X – B _F , i.e, four different problems in all. The functionality of our implementation corresponds to that of the LAPACK routine DGELS. The new implementation is significantly faster and simpler. It outperforms the LAPACK DGELS for all matrix sizes tested. The improvement is usually 50–100% and it is as high as 400%. The four different problems of DGELS are essentially reduced to two, by use of explicit transposition of A. By explicit transposition we avoid computing Householder transformations on vectors with large stride. The QR factorization of block columns of A is performed using a recursive level-3 algorithm. By interleaving updates of B with the factorization of A, we reduce the number of floating point operations performed for the linear least squares problem. By avoiding redundant computations in the update of B we reduce the work needed to compute the minimum norm solution. Finally, we outline fully recursive algorithms for the four problems of DGELS as well as for QR factorization.This revised version was published online in October 2005 with corrections to the Cover Date.     

A Note on the Optimal L 2-Estimate of the Finite Volume Element Method

In this note, the optimal L ²-error estimate of the finite volume element method (FVE) for elliptic boundary value problem is discussed. It is shown that u–u _{h 0}Ch ²|ln h|^1/2f_1,1 and u–u _{h 0}Ch ²f_1,p , p>1, where u is the solution of the variational problem of the second order elliptic partial differential equation, u _h is the solution of the FVE scheme for solving the problem, and f is the given function in the right-hand side of the equation.     

Small into isomorphisms onL ∞ spaces

IfT is an isomorphism ofL (A, ) intoL (B, ) which satisfies the condition T T ^–11+, where (A, ) is a -finite measure space, thenT/T is close to an isometry with an error less than 4.     

On the norm of polynomials of two adjoint projections

Letf(X, Y) be a polynomial of two non-commuting variables and letP be an arbitrary nontrivial projection operator in Hilbert space. The class of all polynomialsf(X, Y) for which f(P, P ^*) depends only onf and P are described. In the case when such a dependence exists the explicit formula is obtained. Some applications to singular integral operators are given.     

An algorithm for best approximate solutions ofAx=b with a smooth strictly convex norm

Summary In this paper, overdetermined systems ofm linear equations inn unknowns are considered. With ^m equipped with a smooth strictly convex norm, ·, an iterative algorithm for finding the best approximate solution of the linear system which minimizes the ·-error is given. The convergence of the algorithm is established and numerical results are presented for the case when · is anl _p norm, 1<p<.Portions of this paper are taken from the author's Ph.D. thesis at Michigan State University     

Characterizations of logA≥logB and normaloid operators via Heinz inequality

In 1951, Heinz showed the following useful norm inequality:If A, B0and XB(H), then AXB ^r X^1–r A ^r XB ^r holds for r [0, 1]. In this paper, we shall show the following two applications of this inequality:Firstly, by using Furuta inequality, we shall show an extension of Cordes inequality. And we shall show a characterization of chaotic order (i.e., logAlogB) by a norm inequality.Secondly, we shall study the condition under which , where is Aluthge transformation ofT. Moreover we shall show a characterization of normaloid operators (i.e.,r(T)=T) via Aluthge transformation.