Roundoff-error analysis of a new class of conjugate-gradient algorithms

We perform the rounding-error analysis of the conjugate-gradient algorithms for the solution of a large system of linear equations Ax=b where Ais an hermitian and positive definite matrix. We propose a new class of conjugate-gradient algorithms and prove that in the spectral norm the relative error of the computed sequence {x_k} (in floating-point arithmetic) depends at worst on ζк^{$\text{3}\text{2}$}, where ζ is the relative computer precision and к is the condition number of A. We show that the residual vectors r_k=Ax_k-b are at worst of order ζк?vA?v ?vx_k?v. We p oint out that with iterative refinement these algorithms are numerically stable. If ζк ² is at most of order unity, then they are also well behaved.     

Bounds for the error of linear systems of equations using the theory of moments

Consider the system, of linear equations Ax = b where A is an n × n real symmetric, positive definite matrix and b is a known vector. Suppose we are given an approximation to x, ξ, and we wish to determine upper and lower bounds for ∥ x − ξ ∥ where ∥ ··· ∥ indicates the euclidean norm. Given the sequence of vectors {r_i}_{i^k = 0}, where r_i = Ar_{i − 1} and r₀ = b − Aξ, it is shown how to construct a sequence of upper and lower bounds for ∥ x − ξ ∥ using the theory of moments.     

A posteriori residual error estimation of a cell-centered finite volume method

We present an a posteriori residual error estimator for the Laplace equation using a cell-centered finite volume method in the plane. For that purpose we associate to the approximated solution a kind of Morley interpolant. The error is then the difference between the exact solution and this Morley interpolant. The residual error estimator is based on the jump of normal and tangential derivatives of the Morley interpolant. The equivalence between the discrete H¹-seminorm of the error and the residual error estimator is proved. The proof of the upper error bound uses the Helmholtz decomposition of the broken gradient of the error and some quasi-orthogonality relations. To cite this article: S. Nicaise, C. R. Acad. Sci. Paris, Ser. I 338 (2004).     

On the uniqueness of overcomplete dictionaries, and a practical way to retrieve them

A full-rank under-determined linear system of equations Ax = b has in general infinitely many possible solutions. In recent years there is a growing interest in the sparsest solution of this equation—the one with the fewest non-zero entries, measured by ∥x∥₀. Such solutions find applications in signal and image processing, where the topic is typically referred to as “sparse representation”. Considering the columns of A as atoms of a dictionary, it is assumed that a given signal b is a linear composition of few such atoms. Recent work established that if the desired solution x is sparse enough, uniqueness of such a result is guaranteed. Also, pursuit algorithms, approximation solvers for the above problem, are guaranteed to succeed in finding this solution.Armed with these recent results, the problem can be reversed, and formed as an implied matrix factorization problem: Given a set of vectors {b_i}, known to emerge from such sparse constructions, Ax_i = b_i, with sufficiently sparse representations x_i, we seek the matrix A. In this paper we present both theoretical and algorithmic studies of this problem. We establish the uniqueness of the dictionary A, depending on the quantity and nature of the set {b_i}, and the sparsity of {x_i}. We also describe a recently developed algorithm, the K-SVD, that practically find the matrix A, in a manner similar to the K-Means algorithm. Finally, we demonstrate this algorithm on several stylized applications in image processing.     

Information Complexity of Multivariate Fredholm Integral Equations in Sobolev Classes

In this paper, the complexity of full solution of Fredholm integral equations of the second kind with data from the Sobolev classW^r₂is studied. The exact order of information complexity is derived. The lower bound is proved using a Gelfand number technique. The upper bound is shown by providing a concrete algorithm of optimal order, based on a specific hyperbolic cross approximation of the kernel function. Numerical experiments are included, comparing the optimal algorithm with the standard Galerkin method.     

On the asymptotic convergence factor of the total step method in interval computation

For a special class of n×n interval matrices A we derive a necessary and sufficient condition for the asymptotic convergence factor α of the total step method x^(m+1)=Ax^(m)+b to be less than the spectral radius ϱ(|A|) of the absolute value |A| of A.     

Linear equations in infinite matrices

The existence and uniqueness of a bounded solution are shown for the linear equation in infinite matrix $\text{A}\text{x}\text{=}\text{b}$ where A is strictly diagonally dominant and b is bounded.     

A computational method for the indefinite quadratic programming problem

We present an algorithm for the quadratic programming problem of determining a local minimum of ?(x)=$\text{1}\text{2}$x^TQx+c^Tx such that A^Tx?b where Q ymmetric matrix which may not be positive definite. Our method combines the active constraint strategy of Murray with the Bunch-Kaufman algorithm for the stable decomposition of a symmetric matrix. Under the active constraint strategy one solves a sequence of equality constrained problems, the equality constraints being chosen from the inequality constraints defining the original problem. The sequence is chosen so that ?(x) continues to decrease and x remains feasible. Each equality constrained subproblem requires the solution of a linear system with the projected Hessian matrix, which is symmetric but not necessarily positive definite. The Bunch-Kaufman algorithm computes a decomposition which facilitates the stable determination of the solution to the linear system. The heart of this paper is a set of algorithms for updating the decomposition as the method progresses through the sequence of equality constrained problems. The algorithm has been implemented in a FORTRAN program, and a numerical example is given.     

Arithmetic behaviour of the sums of three squares

Let n ≠ 4^a(8b + 7) be an integer. We deal with the problem of the solvability of the equation n = x₁² + x₂² + x₃² in integers x₁, x₂, x₃ prime to n. By a theorem of Vila (Arch. Math. 44 (1985), 424–437), the existence of such a solution implies that every central extension of the alternating group A_n, for n ≡ 3 (mod 8), can be realized as a Galois group over Q.     

Optimal approximation of convex curves by functions which are piecewise linear

In this paper an efficient method is presented for solving the problem of approximation of convex curves by functions that are piecewise linear, in such a manner that the maximum absolute value of the approximation error is minimized. The method requires the curves to be convex on the approximation interval only. The boundary values of the approximation function can be either free or specified. The method is based on the property of the optimal solution to be such that each linear segment approximates the curve on its interval optimally while the optimal error is uniformly distributed among the linear segments of the approximation function. Using this method the optimal solution can be determined analytically to the full extent in certain cases, as it was done for functions x² and $\text{x}{}^{\mathrm{<mtext>1</mtext><mtext>2</mtext>}}$. In general, the optimal solution has to be computed numerically following the procedure suggested in the paper. Using this procedure, optimal solutions were computed for functions sin x, tg x, and arc tg x. Optimal solutions to these functions were used in practical applications.     

An Improved Interval Linearization for Solving Nonlinear Problems

Let f(x) denote a system of n nonlinear functions in m variables, m≥n. Recently, a linearization of f(x) in a box x has been suggested in the form L(x)=Ax+b where A is a real n×m matrix and b is an interval n-dimensional vector. Here, an improved linearization L(x,y)=Ax+By+b, x∈x, y∈y is proposed where y is a p-dimensional vector belonging to the interval vector y while A and B are real matrices of appropriate dimensions and b is a real vector. The new linearization can be employed in solving various nonlinear problems: global solution of nonlinear systems, bounding the solution set of underdetermined systems of equations or systems of equalities and inequalities, global optimization. Numerical examples illustrating the superiority of L(x,y)=Ax+By+b over L(x)=Ax+b have been solved for the case where the problem is the global solution of a system of nonlinear equations (n=m).     

Linear equations over commutative rings

A generalized rank (McCoy rank) of a matrix with entries in a commutative ring R with identity is discussed. Some necessary and sufficient conditions for the solvability of the linear equation Ax = b are derived, where x, b are vectors and A is a matrix with entries in either a Noetherian full quotient ring or a zero dimensional ring.     

A new method for computing the matrix exponential operation based on vector valued rational approximations

In this paper a new method for computing the action of the matrix exponential on a vector e^Atb, where A is a complex matrix and t is a positive real number, is proposed. Our approach is based on vector valued rational approximation where the approximants are determined by the denominator polynomials whose coefficients are obtained by solving an inexpensive linear least-squares problem. No matrix multiplications or divisions but matrix-vector products are required in the whole process. A technique of scaling and recurrence enables our method to be more effective when the problem is for fixed A,b and many values of t. We also give a backward error analysis in exact arithmetic for the truncation errors to derive our new algorithm. Preliminary numerical results illustrate that the new algorithm performs well.     

Why Gaussian quadrature in the complex plane?

This paper synthesizes formally orthogonal polynomials, Gaussian quadrature in the complex plane and the bi-conjugate gradient method together with an application. Classical Gaussian quadrature approximates an integral over (a region of) the real line. We present an extension of Gaussian quadrature over an arc in the complex plane, which we call complex Gaussian quadrature. Since there has not been any particular interest in the numerical evaluation of integrals over the long history of complex function theory, complex Gaussian quadrature is in need of motivation. Gaussian quadrature in the complex plane yields approximations of certain sums connected with the bi-conjugate gradient method. The scattering amplitude c ^T A ^–1 b is an example where A is a discretization of a differential–integral operator corresponding to the scattering problem and b and c are given vectors. The usual method to estimate this is to use c ^T x ^(k). A result of Warnick is that this is identically equal to the complex Gaussian quadrature estimate of 1/. Complex Gaussian quadrature thereby replaces this particular inner product in the estimate of the scattering amplitude.     

Structure of A*-Fibrations over One-Dimensional Seminormal Semilocal Domains

In this paper we give a structure theorem for an A*-fibration over a one-dimensional noetherian seminormal semilocal domain and show that, in this situation, any A*-fibration whose spectrum occurs as an affine open subscheme of the spectrum of an A¹-fibration (equivalently, an affine line A¹) is actually A*. The structure theorem provides examples of A*-fibrations over one-dimensional noetherian seminormal semilocal domains whose spectra are not affine open subschemes of any affine line A¹ over the base ring. We also construct examples of nontrivial A*-fibrations over one-dimensional noetherian non-seminormal local domains whose spectra are open subschemes of A¹-fibrations over the base ring.     

L p -Norms, Log-Barriers and Cramer Transform in Optimization

We show that the Laplace approximation of a supremum by L ^p -norms has interesting consequences in optimization. For instance, the logarithmic barrier functions (LBF) of a primal convex problem P and its dual P ^* appear naturally when using this simple approximation technique for the value function g of P or its Legendre–Fenchel conjugate g ^*. In addition, minimizing the LBF of the dual P ^* is just evaluating the Cramer transform of the Laplace approximation of g. Finally, this technique permits to sometimes define an explicit dual problem P ^* in cases when the Legendre–Fenchel conjugate g ^* cannot be derived explicitly from its definition.     

On the rational approximation ofe −x on [0, ∞)

In recent years the problem of uniform approximation ofe ^?x on [0, ∞) by rational functions has found wide interest. In this paper we present a method for the construction of rational approximations toe ^?x that seem to come arbitrarily near to the asymptotically best one. We start with a generalization of the integral form of the Padé approximant by introducing certain real parametersa _j ,b _i ,k and?. The corresponding error function has again an integral representation and is estimated for everyx∈[0,∞) by the Laplace method. This leads to the consideration of a finite number of new error functionsφ _i·j whose maxima are determined by a set of nonlinear equations and, under some restrictions on thea _j ,b _i ,k, and?, are also unique. Variation of these parameters according to a special algorithm and computation of the corresponding maxima of the functionsφ _i·j shows that forn→∞ the order of rational approximation ofe ^?x on [0, ∞) is at least 9.03^?n .     

Uniform approximation by some Hermite interpolating splines

In this paper some upper bound for the error ∥ s-f ∥_∞ is given, where f ε C¹[a,b], but s is a so-called Hermite spline interpolant (HSI) of degree 2q ?1 such that f(x_i) = s(x_i), f′(rmx_i) = s′(x_i), s^(j) (x_i) = 0 (i = 0, 1, …, n; j = 2, 3, …, q ?1; n > 0, q > 0) and the knots x_i are such that a = x₀ < x₁ < … < x_n = b. Necessary and sufficient conditions for the existence of convex HSI are given and upper error bound for approximation of the function fε C¹[a, b] by convex HSI is also given.     

On a viscous Hamilton-Jacobi equation with an unbounded potential term

We present comparison, uniqueness and existence results for unbounded solutions of a viscous Hamilton-Jacobi or eikonal equation. The equation includes an unbounded potential term V(x) subject to a quadratic upper bound. The results are obtained through a tailor-made change of variables in combination with the Hopf-Cole transformation. An integral representation formula for the solution of the Cauchy problem is derived in the case where V(x)=ω²|x|²/2.     

A Dual Approach to Constrained Interpolationfrom a Convex Subset of Hilbert Space

Many interesting and important problems of best approximationare included in (or can be reduced to) one of the followingtype: in a Hilbert spaceX, find the best approximationP_K(x) to anyxXfrom the setKC∩A⁻¹(b),whereCis a closed convex subset ofX,Ais a bounded linearoperator fromXinto a finite-dimensional Hilbert spaceY, andbY. The main point of this paper is to show thatP_K(x)isidenticaltoP_C(x+A*y)—the best approximationto a certain perturbationx+A*yofx—from the convexsetCor from a certain convex extremal subsetC_bofC. Thelatter best approximation is generally much easier to computethan the former. Prior to this, the result had been known onlyin the case of a convex cone or forspecialdata sets associatedwith a closed convex set. In fact, we give anintrinsic characterizationof those pairs of setsCandA⁻¹(b) for which this canalways be done. Finally, in many cases, the best approximationP_C(x+A*y) can be obtained numerically from existingalgorithms or from modifications to existing algorithms. Wegive such an algorithm and prove its convergence